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Abstract 


A study into the problem of determining electromagnetic solutions at high frequencies for 
problems involving complex geometries, large sizes and multiple sources (e.g. antennas) has been 
initiated. Typical applications include the behavior of antennas (and radiators) installed on complex 
conducting structures (e.g. ships, aircrafts, etc..) with strong interactions between antennas, the 
radiation patterns, and electromagnetic signals is of great interest for electromagnetic compatibility 
control. This includes the overall performance evaluation and control of all on-board radiating 
systems, electromagnetic interference, and personnel radiation hazards. 

Electromagnetic computational capability exists at NASA LaRC, and many of the codes 
developed are based on the Moment Method (MM). However, the MM is computationally intensive, 
and this places a limit on the size of objects and structures that can be modeled. Here, two approaches 
are proposed: (i) a current-based hybrid scheme that combines the MM with Physical optics, and (ii) 
an Alternating Direction Implicit-Finite Difference Time Domain (ADI-FDTD) method. The essence 
of a hybrid technique is to split the overall scattering surface(s) into two regions: (a) a MM zone 
(MMZ) which can be used over any part of the given geometry, but is most essential over irregular 
and “non-smooth” geometries, and (b) a PO sub-region (POSR). Currents induced on the scattering 
and reflecting surfaces can then be computed in two ways depending on whether the region belonged 
to the MMZ or was part of the POSR. For the MMZ, the current calculations proceed in terms of 
basis functions with undetermined coefficients (as in the usual MM method), and the answer obtained 
by solving a system of linear equations. Over the POSR, conduction is obtained as a superposition of 
two contributions: (i) currents due to the incident magnetic field, and (ii) currents produced by the 
mutual induction from conduction within the MMZ. This effectively leads to a reduction in the size 
of linear equations from N to N - Npo with N being the total number of segments for the entire 
surface and Npo the number of segments over the POSR. The scheme would be appropriate for 
relatively large, flat surfaces, and at high frequencies. 

The ADI-FDTD scheme provides for both transient and steady state analyses. The restrictive 
Courant-Friedrich-Levy (CFL) condition on the time-step is removed, and so large time steps can be 
chosen even though the spatial grids are small. This report includes the problem definition, a detailed 
discussion of both the numerical techniques, and numerical implementations for simple surface 
geometries. Numerical solutions have been derived for a few simple situations. 
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I. Introduction and Problem Statement : 

Research interest at NASA in the general area of electromagnetics arises from a broad range 
of potential applications involving detection and remote sensing, system stability with regard to 
electromagnetic disturbances/pulses, stochastic variability, and scattering effects in radar systems. 
Many of these applications rely on the use of antennas that may be mounted on or placed in close 
proximity to suitable platforms (for example, airplanes). It therefore, becomes important to be able to 
accurately analyze and predict antenna performance characteristics in the presence of such platforms 
and large objects involving complex geometries. 

Electromagnetic solutions (including parameters such as impedance and radiation patterns) 
can be obtained by either the moment method [1,2], or finite-element analysis [3]. Of these, the 
moment method (MM) is more accurate, but computationally intensive as it leads to a dense matrix 
that has to solved. The finite element method (FEM) is more attractive since it results in sparse, 
banded matrices which can stored and solved more easily. However, the FEM for unbounded 
problems does not include the Sommerfield radiation condition, and requires the use of adequate 
absorbing boundary conditions to treat the radiation conditions that must be imposed at the 
boundaries of discretized region. The MM scheme is accurate but slow and computationally very 
intensive, especially as the physical dimensions approach lengths on the order of ten wavelengths. 
Moreover, MM is applicable to complex conducting or penetrable bodies, and complicated 
geometries. The MM requires the solution of “N” equation, where “N” the number of basis functions 
chosen depend on the number of patch-elements (or segments) into which the overall physical system 
and structure is decomposed. These “N” equations have to be solved, which typically requires storage 
on the order of N 2 and computation times on the order of N k , with 2 < k < 3 [4]. In employing basis 
functions, the segment length has to be taken to be smaller than A/10, with A being the operating 
wavelength. Obviously, for high-frequency situations, A reduces dramatically, leading to increases in 
both the storage and computational requirements which scale as f 4 and f 2k . Consequently, from a 
practical standpoint, the use for the MM technique is restricted to applications involving low 
frequencies and/or medium sized scattering bodies. At high frequencies, or for large sized bodies, 
the MM becomes computationally intractable. 

The problem of determining electromagnetic solutions at high frequencies can be simplified 
by incorporating ideas of “optics”. An important premise for the optics-based techniques is that at 
sufficiently high frequencies the propagation, scattering and diffraction of electromagnetic signals 
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exhibit highly localized behavior. Such a localized behavior can qualitatively be demonstrated by 
considering the radiation integral over the spatial current distribution induced on a radiating object by 
the primary excitation. At high frequency, the integrand oscillates rapidly to produce destructive 
interference. The dominant contribution arises only from waves that emanate from the local 
neighborhood of critical areas that are points of reflection, transmission and diffraction [5]. The 
connection between ray optics and Maxwell’s equations has also been demonstrated by deriving the 
laws of geometric optics from Helmholtz’s scalar wave equation [6,7]. In such a simple ray treatment 
based solely on incident, reflected and transmitted beams, no rays can exist in the “shadow 
boundaries” beyond the “line-of-sight” visibility regions. 

An ideal option would be the use of hybrid schemes that preserve the accuracy without 
computational complexity, as suggested by Thiele et al. [8], Burnside and co-workers [9], and 
reviewed more recently [10,1 1 ]. In general, the hybrid scheme can be formulated as a field-based 
(FB) or current-based (CB) analysis. In the field-based version [12,13], solutions for the fields 
associated with edge or surface diffraction, based on the geometrical theory of diffraction (GTD) or 
geometrical optics (GO), could be used as the starting point. These could serve as an Ansatz to the 
MM formulation and be a representation for parts of a scatterer not conforming to a canonic 
geometry. The FB approach is attractive for problems where the radiator is in proximity to large 
convex surfaces for which the canonical GTD solutions are known. These approaches are also known 
as “ray-based” schemes provide very fast results, but their application is restricted to problems where 
asymptotic high-frequency solutions are available. In these GTD-MM schemes, the form of the 
currents induced by straight-edge diffraction or in the shadow regions were derived from GTD. In 
CB formulations, the analysis proceeds from Ansatz solutions for the currents obtained from optics- 
theory [14], such as physical optics (PO) or the Fock theory. A Galerkin representation of the 
unknown currents is used in regions of the scatterer where Ansatz solutions do not exist, such as 
surfaces with materials and/or geometric irregularities. The CB approach can be more advantageous 
since the MoM is based on currents as well. Hence, a continuous current flow can be modeled on the 
whole surface of the scattering body. Also, since CB-based methods are inherently capable of 
modeling arbitrary-shaped geometries. A number of CB hybrid methods have been discussed in the 
literature [4, 15-17]. 

Here, a hybrid technique that combines the MM with aspects of PO has been used. It requires 
less computational work than a pure MM scheme. Such as approach would be suitable for high- 
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frequency analysis or for situations in which the physical dimensions of scattering objects was large. 
The essence of a hybrid technique would be to invoke a splitting of the overall scattering surface(s) 
into two regions: (i) a MM zone (MMZ) which can be used over any part of the given geometry, but 
is most essential over irregular and “non-smooth” geometries, and (ii) a PO sub-region (POSR). 
Currents induced on the scattering and reflecting surfaces could then be computed in two ways 
depending on whether the region belonged to the MMZ or was part of the POSR. For the MMZ, the 
current calculations would proceed in terms of the basis functions with coefficients determined 
through the solution of a linear system of equations (usual MM method). Over the POSR, a 
superposition of basis functions with known coefficients could be used. This could effectively lead 
to a reduction in the size of linear equations from N to N - Npo with N being the total number of 
segments for the entire surface and Npo the number of segments over the POSR. However, the 
approximations inherent in such as hybrid approach are listed below. For self-consistency and 
accuracy, both the self- and mutual-coupling between the MM and PO regions have to be considered. 

Assumptions Inherentjn the Hybrid, MM-PO Approach: 

(i) PO approach in the illuminated region, and for a perfect electric or magnetic conductor. 

(ii) Secondary radiation due to many-body inter-coupling, creeping wave effects and edge 
diffraction neglected. Such effects could only be accounted for by combining the MM with 
other techniques such as GTD, Fock theory etc.. 

(iii) PO applied only over regions of “smooth” geometry. All other irregular or “non-smooth” 
geometric regions assigned to the MM domain. 

(iv) Hybrid approach only valid for high frequencies or when dimensions of physical object much 
larger than operating wavelength. 

(v) Applicable only for linear systems i.e. media for which the permeability (u), permittivity (e) 
and conductivity (a) are constant and not functions of the field excitation. 

(vi) Accuracy limited by the finiteness of the basis functions - since an exact MM solution is only 
available in theory for an infinite basis set. 

Advantages of_ the ADI-FDTD Simulation Approach: 

(i) The ADI-FDTD simulation scheme allows for both a complete time-domain analysis and 
the steady-state result. The latter can result by running the simulations up to longer times. 
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The transient behavior cannot be obtained from the Moment Method or PO-MoM 
schemes. 

(ii) The ADI-FDTD allows for any arbitrary geometry including sharp comers and complex 
features. 

(iii) It is significantly faster than the traditional FDTD scheme since very large time steps can 
be used without the need to be restricted to the CFL condition. 

(iv) High frequency analyses for which the spatial discretization has to be small, can be 
carried out relatively easily as the time-step can remain large. Similarly, the technique is 
of advantage for large objects for which many grid points may be required. 

(v) Non-uniform, multi-grids can be employed for better resolution. 

(vi) Absorbing boundary conditions can be implemented for accurate analyses. 

(vii) Finally, dielectrics or perfect conductors or both can be studied with ADI-FDTD. Thus, 
if becomes easy to look at problem that involve coatings, and multiple layers. 

II. Theoretical Background of the Moment Method : 

A time-harmonic electric field E inc (r) incident on a conducting structure induces currents 
J(r’) on the surface, which in turn, radiate a secondary (or scattered) field E ,ca (r) given by: 

E* c “ (r) = - j w A(J, r) - V4>(p,r) , (1) 

where A(J,r) the magnetic vector potential is: 

A(J,r) = [p/(4n)j J [{J(r’) exp[-jkR]}/R ] dV' , (2a) 

and a>(p,r) the electric scalar potential is expressed as: 

<t(p,r) = [l/(4ne)]/[p(r’)exp[-jkR]}/R] dV' , (2b) 

with R being the distance from the source point (r’) to the observation point (r), i.e. R = |r - r*| . 
Since the tangential component of the total electric field must equal zero at the boundary of a perfect 
conductor, one gets: E‘°‘ (r)U = E sca (r)^ + E lnc (r)|tan = 0 E ,nc (r)Un = -E sc<, (r)| tan . Using eqns. 
(1) and (2a,2b), this yields the following electric field integral equation (EFIE) for the unknown 
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current density J(r’) : 

E toc (r)| tan = Ljoju/(4n)J [ { J(r’)exp[-jkR] }/R]dV -l/(4njo)e)] V{JV.J (r») exp[-jkR]/R}dV]U , (3) 

Where the continuity equation: p(r’) = -1 /(jco) V.J (r’) has been used to relate the charge density 
p(r’) to the current density J (r’). 

In this EFIE formulation, the current density J (r’) is the unknown to be solved, given a 
complete specification of the incident electric field E lnc (r). The EFIE, which a linear integro- 
differential equation has the following general form: Lf(x) = g(x), where g(x) denotes the incident E- 
field excitation and f(x) the unknown currents. The moment method (MM) consists of converting the 
EFIE into a matrix equation that can be solved for the unknown currents [1]. This is done by 
expanding the unknown f(x) in terms of a finite number of sub-domain basis function f B i(x) as: f(x) = 
i£ n a; f B i(x) . Using these basis functions in the EFIE, and taking an inner product with of both sides 
with respect to testing functions fn(x) converts the integral equation into a system of N equations in 
the N-unknowns a* . Thus, for example, 

i£ N ai < f T j(x), Lf Bi (x) > = < f Tj (x) , g(x) > , j = 1, 2, .... N , (4) 

i.e. [Lji] [ai] = [gj ] . 

By determining the unknown coefficients aj , the desired solution f(x) to the EFIE is obtained. 
Usually, the test functions are taken to be the same as the basis functions [i.e. fn(x) = f B i(x)] to ensure 
that the matrix [L] is nonsingular. 

II.A Simple Planar Surfaces for the Moment Method : 

In general, the electromagnetic scattering problems of interest present the following four types 
of structures [18]: (i) Surfaces - planar or curved, (ii) wires as in current driven wire antennas, (iii) 
wire-wire junctions, and (iv) surface-wire junctions. The basis functions B„(r) for surfaces were first 
proposed by Rao et al. [19] and describe current flow on a surface between triangular patches. One 
basis function is associated with the n* edge and is defined to be: 
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B„(r) = {Ln /(2A n± )} p„+ , for r e T n± 

= 0 , elsewhere , (5) 

where Ln is the common edge between two triangles T n + and T n „ , having areas A n± , respectively. 
The global position r and the local position vectors p„+ are shown in Fig. 1 . The reference direction 
of current flow is from T n + to T n . and the total charge associated with each pair of triangles is zero 
[ 18 ]. 

Basis functions for wires, wire-surface junctions and wire- wire junctions have similarly been 
determined in the literature, and can be found in a review report by Veihl and Mittra [1 8]. However, 
the details are not given here in this interim report. The omission arises from a desire to keep the 
structures simple during this initial developmental stage, and because the MM-PO formulation has 
only been applied thus far to surfaces. Based on the reports in the literature (e.g. [18]), In the near 
the calculations can easily be extended to include other physical structures (such as wires and wire- 
surface junctions) in the near future. 

Using the surface basis functions, as a concrete yet simple demonstrative example, the current 
J(R) can be expressed in terms of the bases Bi(r) as: J(R) = i£ N Ii(r) Bi(r) . It must be mentioned, 
that here again for simplicity, only a single surface has been considered. For multiple surfaces, the 
above single sum will transform over to a double summation. Consequently, the scattered electric 
field E“*(r) can be expressed in terms of the vector and scalar potentials and the bases functions as: 

E ,C V)= iL N IiEi(r) , (6a) 

where, Ei(r) = - j q Ai(r) - V<l>i(r) , (6b) 

A|(r) = pf/s- Bi(r’) G(r,r’) dS’ , (6c) 

and, <J>i(r) = (1/e) JJs> b,(r’)G(r,r’)dS’ , (6d) 

where as given in [18], bi(r’) = -[l/(j a)] V. B|(r’) . For the surface basis functions of eqn. (5), 
bi(r’) = - { ± Ln /(jcoA n +)} for r’ e T n + , and bj(r’) = 0 elsewhere. In eqns. (6c) and (6d) above, 
G(r,r’) is the free space Green’s function given by: 

G(r,r’) = exp[-jkR]/(4nR) , (6e) 

with R = |r-r’| . Substitution of the Green’s function and the surface basis function Bi(r’) from 
eqn. (5) into eqns. (6) leads to the following expressions for the vector and scalar potentials: 
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Fig. 1 . Definition of the surface basis functions B„(r) in terms of the global co-ordinates associated 
with the origin O. The triangles T n+ and T n ._ on either side of the common boundary (i.e. the n* 
edge) having length U are shown [after reference 18]. 
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A,(r) - {ia/(4n)}[J7s'+ [Li/(2A i+ R)] exp[-jkR]p 1+ . dS’+ f JV- [L i /(2A i .R)]exp[-jkR]p i _dS’] , (7a) 


and, <X>i(r) = {-l/(4njcoe)}[JJV + [L,/(A i+ R)] exp[-jkR] dS’- J/s-. [L i /(A i .R)]exp[-jkR]dS’] . (7b) 

Since E lnc = -E ,ca , the governing MM equation becomes: <E tac ,B n > = i£ n Ii < B n ,(jcoA ( + V<J>i ) > . 
The use of (5) for the bases vectors and the one-point approximation for the surface integral at the 
centroid, leads to: 

< E lnc . B n > = (Ln/2) [ E inc (r cn+ ) . p cn+ + E ,nc (r cn .) . p cn . ] , (8a) 

^ B n .(JcoAi ) > — (jcoLn/2) [Ai(r cn +) • Pcn+ A|(r cn .) . Pen- ] > (8b) 
and , < B n . > < (V . B n ) - > — Ln [^i (l"cn-) - (l*cn+)] • (8c) 

This can be cast in the usual matrix notation as: [V] = [Z] [I] , where the elements of the matrix [V] 
are: V„ = < E lnc . B n > , and the elements of [Z] are: - <B„ .(jcoAi ) > + < B„ .V<t>j > . The 

evaluation of V n is straight-forward based on eqn. (8a). Hence, from a numerical implementation 
standpoint, the main task in the solution of the system of linear equations is the evaluation of the 
matrix elements Zn, . The computation for Z™ are discussed next. 

II.B Evaluation of the Matrix Elements : 

Inspection of eqns. (8b) and (8c), reveals that the matrix elements Z„j contain two types of 

terms: Ai(r c „).p cn and $i(r c „) . The A|(r cn ).Pcn terms require surface integrations Ii of the form: 

I, ~ [Li/(2A0] { J/s- (1/R) exp[-jkR] p, dS’ } , (9a) 

where R = | r cn - r’ | . Here, r’ is a source point in the i* triangle, pi the position vector from the 
vertex of the i A triangle to the source point, and r cn the centroid of the n* 11 triangle. Similarly, all of 
the <l>i (r c „) terms require surface integrations I 2 of the form: 

h ~ [Li/A,] { J/s- (1/R) exp[-jkR] dS’ } . (9b) 
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If the 11 th and i** 1 triangles are different (i.e. for n * i), the value of R is always non-zero, and no 
singularity results. In such cases, the following methods can all be used for evaluating the triangular 
double integrals: (i) A one-point approximation by sampling at the centroid of the triangle. This is 
the least accurate method, but could suffice if R » Rcmax , with Rc max being the maximum distance 
between the centroids of all the patches, (ii) Use of a seven-point numerical integration valid for 
triangular sources as given by Stroud [20]. (iii) A numerical integration technique based on the 
normalized area co-ordinates [21] as discussed by Rao et al. [19]. However, if i - n, then R can 
approach zero and the integrands become singular. For such situations, the singularity has to be 
removed by decomposing the integrals into an analytical and a numerical part as outlined in Ref. 1 8. 
Thus: 

I, ~ { fJV (l/R) exp[-jkR] pi dS’ } = JJ s -(l/R){exp[-jkR] -1 }p,dS’ + JJ S '(pi/R) dS’ , (10a) 
and, I 2 ~ { JJ S - (l/R) exp[-jkR] dS’ } = J/ s <l/R){exp[-jkR] -1 }dS’ + JJV(1/R) dS’ , (10b) 

i.e. Ii = In + 1] 2 , and I 2 = I 2 i + 122 • Of the four integral terms on the right sides of eqns. (10a) and 
(10b), In and I 2 i can both be integrated numerically, since LtR_>o [{exp[-jkR] — 1 }/R] 4 -jk, and 
hence, the integrand is finite and non-divergent. The integrals Ii 2 and I 22 can be evaluated 
analytically[18, 21,22], and the resulting expressions can best be understood in terms of Fig. 2 
showing a planar triangle. In Fig. 2, the triangle is defined by a local reference frame (u,v,w) and 
assumed to lie in the w = 0 plane. Each side of length Lj (i= 1,2,3) is shown to be opposite node ni 
with co-ordinates (0,0,0), (L 3 ,0,0) and (u 3 ,v 3 ,0). Also, as shown in Fig. 2, tj + and tf are distances 
from the observation point (uo,vo,0) to the end points of each side Li . From a geometric analysis, one 
obtains : 

si' = - [(L 3 -uo)(L 3 -u 3 ) + v 0 v 3 ]/Li , (11a) 

s, + = [(u 3 - u 0 ) (u 3 - L 3 ) + v 3 (v 3 - v 0 )] /Li , (lib) 
s 2 * = - [u 3 (u 3 - Uo) + v 3 (v 3 - v 0 )]] /L 2 , (lib) 
s 2 + = [u 0 u 3 + v 0 v 3 ] /L 2 , (lid) 

s 3 ‘ = -uo, ands 3 + = L 3 -u 0 , (1 le) tf = [(L 3 - uo) 2 + (v 3 - v 0 ) 2 ] 05 , (Ilf) 

t, + = [(u 3 - Uo) 2 + (v 3 - v 0 ) 2 ] 05 , (llg) 
t 2 + = [uo 2 + v 0 2 ] 0 5 , (llh) 
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(Uq.Vq.O) =p 


Fig. 2. Triangular geometry for the integral evaluations [after reference 22]. (a) Plane triangle T in 
its local frame (u,v,w = 0). The length of the i* side 87 { is denoted by Lj . (b) The tangent s A and 

normal m A unit vectors on the triangle contour, (c) The distances t, and Sj defined in the text from the 
local point at (uo , Vo , 0). 



t3 + - tf ; t2 -ti + ; t3 - 12 + (Hi) 

t,° = Vo [(U3 - U) + V3 (L 3 - uo)] / L, , (1 lj) 

t 2 ° = [ uo u 3 - v 0 113)] / L 2 , and t 3 ° = v 0 . (Ilk) 


In terms of the geometric parameters defined in eqn. (1 la-1 lk), the integrals I 12 and I 22 are [18,22] : 

I 12 = J/s- ((P^) dS’ = I + p cn I 22 , (12a) 
where , I = 0.5 i=1 £ 3 m, {(t, 0 ) 2 Ln [{tj + + Si + }/{tf + sf}] + t*V- tf sf } , (12b) 
and, I 22 = JJs-O/R) dS’ = i=i^ 3 tj° Ln [ {t, + + Si + }/ {tf + sf}] . (12c) 

In eqn. (12b), m, are the outward unit normals from the lines opposite the i th node as shown in 
Fig. 2. 

II.C Example of Incident Plane Wave for the Moment Method : 

Incident excitation on the planar surface used for the MM above can have a variety of form 
depending on the source. These could include plane waves for scattering problems, delta-gap and 
magnetic frills for linear antenna sources, magnetic currents excitations for slot and aperture 
antennas. Here, for simplicity once again, an incident plane wave excitation has been chosen as a 

demonstrative example. Thus: E lnc (r) = [E e a e + E 0 ] exp[-jk.r] , with ae and a$ being unit 
vectors. In terms of Cartesian co-ordinates, these vectors can be expressed as: 

a e = cosQ cos0 a* + cos0 sinp a y - sin0 a z ; a$ = - sincp a x + cos$ a y - sin© a z , (13a) 

and , k = |k| [sin0 cos$ a* + sin0 sincp a y +cos0a z ] . (13b) 


III. Aspects of Physical Optics and the Hybrid MM-PO Scheme : 

As already mentioned in the introduction section, the behavior of antennas (and radiators) 
installed on complex conducting structures (e.g. ships, aircrafts, etc..) are modified with respect to 
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free-space operation due to the interactions between each antenna and the supporting structure. The 
development of numerical methods for estimating the interactions between on-board antennas and 
arbitrary-shaped, on-board, sub-structures is of great interest for electromagnetic compatibility 
control. This includes the overall performance evaluation and control of all on-board radiating 
systems, electromagnetic interference, and personnel radiation hazards. 

Since the MM is computationally intensive and hence, limits the size of objects that can be 
modeled, hybrid techniques that combine the MM with other approximate solutions have begun to 
receive attention [10,1 1,23]. These can be classified as either ray-based [8,9,24] or current-based 
[14-16,25,26] techniques. Of these, the ray-based techniques provide a considerable speed 
advantage, but are difficult to implement for arbitrary geometries and complex objects. By contrast, 
the current-based methods, in which one attempts to determine equivalent surface currents that 
represent the actual objects, are capable of modeling irregular and complex geometries. 

Here, the focus is on time-harmonic problems in which a large perfectly electric conducting 
(PEC) body (e.g. aircraft) has a set of on-board antennas on its surface. This basic problem can be 
replaced by an equivalent problem in which the PEC surface is replaced by equivalent currents 
radiating in space. In the simplest sense, the problem reduces to obtaining the surface currents J 
induced by the excitation fields (E lnc , H lnc ). Since current based MM scheme is computationally 
intensive, a possible solution being attempted here is the use of an alternate hybrid method that 
combines physical optics (PO) with the MM approach. This could be done by dividing the overall 
surface into an MM zone (MMZ) which can be used over any part of the given geometry, but is most 
essential over irregular and “non-smooth” geometries, and the remaining PO sub-region (POSR). As 
discussed in the preceding section, the current J over the MMZ would be given as: J = i£ N I| B|(r) . 
Based on the coefficients I| (obtained from a MM solution over MMZ), the scattered E-field could be 
evaluated using eqn. (6). Such a simple solution, however, would include only the coupling between 
the source radiation and the MMZ, and would not contain any interactions with the POSR. 

The computational scheme being proposed for use here would couple the MM and the POSR 
to include mutual interactions. It is similar to the hybrid methods that have been proposed by Hodges 
and Rahmat-Samii [ 1 6] and Obelleiro et al. [27], Let Ji and J 2 represent the surface currents across 
the MMZ and the POSR sub-domains, respectively. Quite generally, the electric-field integral 
equation (EFIE) can be applied to the MMZ with J 2 taken to be an effective source current, while the 
magnetic-field integral equation (MFIE) applied to POSR with J 2 as the source current. This leads to 
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the following system of coupled equations [16] : 


L e [Ji(r)] + L e [J 2 (r)] = - E inc , tan (r) , for r e MMZ (14a) 
and, J 2 (r) = 2 n A x H inc (r) + L h [Ji(r>] + L h [J 2 (r>] , for r e POSR (14b) 

where n A is the unit surface normal, and the operators L e [J(r)] and Lh[J(r)] are defined as : 

L e [J(r)] = -jcoA(r) - V<t(r) = - jco pj J s - J(r’)G(r,r’)dS’ + {l/(ja>e)}V( J S ’ V . J(r’) G(r,r’) dS\ 

(14c) 

and , L„[J(r)] = 2 n A x JJV J(r’) x [VG(r,r’)] dS’ , (14d) 

with G(r, r’) = {exp[-j k |r- r’|] }/ {4 n |r- r’|} . (14e) 

Eqn. (14b) is complete, yet quite complicated since the current J 2 (r) has three components - (i) the 

n A x H ,nc (r) term originates includes the surface currents induced by the incident magnetic field (due 

to the true sources), (ii) a L h [Ji(r)] term that arises from the mutual interaction effect of currents 
Jj(r) flowing in the MMZ. Such MMZ currents give rise to a scattered electric field, which in turn 
can induce a current component in the POSR, and (iii) the Lh[J 2 (r)] term which denotes self- 
induction contributions. These are the currents induced in the POSR due to currents already flowing 
within this region. 

For a PO implementation, the following approximations are made for simplicity : 

(i) The first term of eqn. (14b) is taken to equal: 2 n A x H ,nc (r) in the “lit” or “directly 

illuminated” parts of the POSR, and is set to zero over the “shadow” sub-sections. The PO approach 
uses the following approximation for current J on a conducting surface with an incident magnetic 

field Hj and unit normal n A : J = 2 y n A x Hi where y is a coefficient that accounts for shadow 
effects. For observation points r lying in the shadow region, y = 0 and y - 1 otherwise. Thus, it is 
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inherently assumed that edge-diffraction and multiple-scattering effects are negligible. Clearly, this 
would be an excellent approximation for surfaces that were large and relatively “flat” . 

(ii) Next, the PO implementation assumes that self-induced currents can be neglected in 
comparison to those arising from the direct incident excitation and mutual-coupling from the MMZ. 
Consequently, the L h [J 2 (r)] term is set to zero, and one effectively obtains the following 
approximation: 


J 2 (r)~ 2 n A x H inc (r) + L„[Ji(r)] , for r e “lit portion of POSR” (14f) 

J 2 (r) ~ Lh[Ji(r)] , for r e “dark/shadow portion of POSR” . (14g) 

The resulting coupled system of equations to be solved in the hybrid MM-PO scheme are thus given 
by eqn. ( 1 4a), eqn. ( 1 4c), and eqn. ( 1 4f-g). Since the incident source field is always specified in any 
problem, the only unknowns involve the currents J j(r) over the MMZ. These can be solved using the 
MM matrix procedure outlined in the previous section. Consequently, the overall solution process 
does not become any more complicated than that for the MM. Thus, by choosing as small an MMZ 
as possible (i.e., including only the irregular or complex geometric portions, or the near- field 
regions), the computational requirement can be greatly reduced. 

III.A Matrix Elements for the Hybrid MM-PO Scheme : 

As with the previous discussion on the MM implementation, the relevant matrix equations for 
the hybrid MM-PO scheme are given here for a simple surface geometry, neglecting any wires, wire- 
surface and wire-wire interactions. In general, the MM could be applied to all wire sections, and 
hence, all wires and wire-sections included in the MMZ. Using the surface basis functions B„(r) as 
given in eqn. (5), the appropriate matrix elements for the hybrid scheme can be obtained. The details 
are given next. 

Using Ji(R) = i£ N Ii(r) Bi(r), eqn. (14a) can be re-written with the aid of eqn. (14c) as: 

-juu/JsMM {i£ N If(r )B|(r’)}G(r,r’)dS’ + {l/(jwe)}Vjf S MM V’ • {,£ N I,(r’)B,(r’)}G(r,r’) dS’ + 
+ L e [J 2 (r)] = - E ln£ |, an (r) , for r e MMZ (15a) 
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where the surface integrals are over the MMZ. Hence, the inner product <E inc |t ail (r).B n (r)> can be 
expressed as the sum of three terms i.e., <E inc |i an (r).B n (r)> = Ti n + T 2n + T 3 n , where : 

<E lnC | tan (r).B n (r)> = (L„/2) [ E ,nc (r cn+ ) . p cn+ + E ,nc (r cn .) . p cn . ] , (from 8a) 

Ti„ = Vi < B„(r ). JJsmm {.Z N Ii(r )B,(r’)}G(r,r’)dS’> , (15b) 

T 2n = - {l/(jcos)} < B n (r ).VJ/smm V* . {,£ N I,(r’)B,(r’) }G(r,r’) dS’> , (15c) 

T 3n = -<B n (r).L e [J 2 (r)]> . (15d) 

Expressions for the three terms Ti n , T 2n and T 3n can now be worked out. Thus: 

Tin — iL If [Aj(r cn +) . pcn+ A,(r cn -) . p C n- ] } > (15c) 

T 2n = .E N I. Ln [®i (r c „.) - 4>i (r cn+ )] , (1 5f) 

T 3 „ = - < B n (r ). {-jcou/Js-J 2 (r’)G(r,r’)dS , + {l/(jcoG)}VjJs-V’.J 2 (r’)G(r,r’)dS’}> , (15g) 

where the surface integration for T 3n is over the POSR. The terms Ti„ and T 2n are the same as 
obtained before in eqn. (8b) and (8c). Using eqns. (14f) and (14g), the T 3n expression can be split 

into two parts T 3n i and T 3n2 with T 3n i containing the 2 n A x H ,nc (r) factor, and T 3n2 the Lh[Ji(r)] 

term. In turn, T 3n i will have two terms due to the two integrals inherent in eqn. (15g). Expressing 
T 3n i = T 3n ] A + T 3n iB to denote the two terms, one has : 


T 3 ni = T 3n , A + T 3n , B - - < B n (r ). {-jcoujjspo {2 n A x H ,nc (r) 1 G(r,r’)dS’> - 
- <B„(r). {l/(jcoe)}V//sp 0 V • {2n A x H inc (r’)}G(r,r’) dS’} > . (15h) 

Hence, 

T 3 nlA — {j^ Lti k'(4n) ( [ pn+(r C n+ )• AV n i+ Pn-( r cn- )• AV n ]„] , (15i) 
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where AV nl± = {J/spo {n A x H inc (r’) } exp(-jk|r cn± - r’ |)/(| r cn± -r’|) dS’ . (15j) 


Similarly, for T 3 n iB one gets : 


T 3 nlB — {* Ln /(j4nC0e)} [ Pn+(l"cn+ )• AVn 2 + + pn— (rcn— )• A Vj^.. ] , (15k) 

where AVn 2 ± = V{/ Jspo V’ . {n A x H inc (r’)} exp(-jk|rc„ ± - r’ |)/(| r cn± -r’|) dS’ . (1 51) 


Thus, both the T 3 n i A and T 3 n iB terms do not depend on the unknown coefficients Ij . Next, from 
eqn. (15g), the T 3n 2 term has the form : 


T 3 n2 = T 3 n 2 A + T 3 n2B = - < B n (r ). {-jco u/Jspo { U[Ji(r)]} G(r,r’)dS’> - 
- <B n (r).{l/(jcoG)}VjJ S PoV’.{L h [J 1 (r)]}G(r,r’)dS’}> . (15m) 

Using the expansion Ji(r) = i£ N Ij Bj(r) , the T 3 n 2 A and T 3n 2 B factors can each be cast as an 
appropriate summation over N weighted by the unknown coefficients I| . Thus: 


T 3 n 2 A = iE N (I,ja)pLn) {[p„ + .(f S po {n A x f/i B|(r”) x V G(r ,r”) dS”} G(r,r’)dS’] + 

+ [p„-.JJspo {n A x J/i Bi(r”) x V G(r ,r”) dS”} G(r,r’)dS’] , (15n) 

where the innermost integral is over the two triangles Tj+ in MMZ, and the outer integral is over the 
entire POSR. Assuming that the MMZ and POSR are well separated so that the one-point formula 
can be used as an approximation for the surface integral, one gets from eqn. (15n) : 


T 3 n2A ll {(Ilj^P-Ln Lj)/(32n )} [Pn+* C|+ "t" Pn+» C|_ + Pn— • Ci+ "I" Pn— • ®1— ] » (15o) 

with c l± = / Jspo {n A x p 1± (r d± )xV’ [exp(-jk|r’- r cl± |) /|r’- r cl± |] {exp(-jk|rc„- r’|) /|r cn - r’| }dS’ , (15p) 
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Similarly, one obtains for T 3 n 2 B : 


T 3 n 2 B — i£ N {“(Ii Ln Lj)/(32n 2 jcos) } [pn+« d|+ + Pn+» d|_ + p n - • d|+ + p n - • d|_ ] , (15q) 

d|+ = V /J S poV’.{n A X Pi+(r C i+)xV’ [exp(-jk|r- r cl± |) /|r- r d± |]{exp(-jk|r cn - r’|) /|r cn - r’| }dS’ , (15r) 

Putting together the Ti n , T 2 n , T 3n iA , T 3n iB , T 3n 2 A and T 3 n 2 B terms a composite matrix equation of the 
form [V + 5V] = [Z + 5Z] . [I] , i.e. [V* ] = [Z* ] . [I] results. Using the following notation for 
compactness: 


P n = (Ln/2) [ E lnc (r cn +) . p cn+ + E inc (rcn-) . p c „- ] , (16a) 

Qn = {j<^ Ln y/(4n)}[ pn+( r cn+ )» AV n l+ + Pn— (*"cn— )• AV n |__ ] , (16b) 

Pn ~ {" Ln/(j4nwe)}[ Pn+(**cn+)« AV n 2 + + Pn— (r c n— )• A Vn 2 ~ ] , (16c) 

S»i ~ {(jwLn/2) [Ai(r cn+ ) • Pcn+ "t" A|(r cn .) . Pen- ] } » (16d) 

Tni — Ln [P] (r cn .) - < l ) i(ren+)] > (16e) 

LJni {(j^l-lLn Li)/(32n )} [Pn+* C|+ "t Pn+. + Pn— • ^1+ Pn — • ^i— ] , (16f) 

V ni = {-(Ln L;)/(32n 2 jcoe)} [p n +. d i+ + p n +. di_ + p„_. d 1+ + p n _. di_ ] , (16g) 

one has : 

[V + 5 V]n = [V* ]„ = [ Pn + Qn +Rn ] , ( 1 6h) 
and , [Z+ 5Z]ni = [Z ]„i = [SW + T ni + U ni + V„] . (16i) 

Solution of the modified matrix [V* ] = [Z* ] [I] , then yields the unknown coefficients [I] as 
before, but with much less computational burden. 

IV. The Finite-Difference, Time-Domain (FDTD) Method 

The finite-difference, time-domain (FDTD) technique, was first proposed by Yee [28] as a 
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numerical method for the direct solution of Maxwell’s time-dependent curl equations. Its popularity 
continues to grow as the computational costs decline. The Maxwell’s curl equations are [29]: 


dt £ £ 

^ =-— VxE-^-H 


(17) 


dt p p 

where o is electrical conductance; p' the magnetic conductance is zero for the present problem. 

In this method, the electric (E) and magnetic (H) fields are discretized across separate nodes 
that are off by half a grid spacing. The time instants at which the E and H fields are updated are also 
off by half a time step. The discretization of the Maxwell curl equations results in six explicit finite 
difference equations [30]. For example, for a two-dimensional (2D) transverse magnetic case, the z 
component of the electric field E z and y component of the magnetic field H y are given as : 

! , a ',j At 1 + ^— A* 4y 


2e i 


hj 


2f, 


*>j 


(18a) 


H n ; V2 (ij ) = h r l/2 (/, j ) - — [ • g " - 1 - ^ — } - ) ] , 

P Ay 


(18b) 


and, H n x +U2 (/, j) = HT xn (/, j) - . 

p Ax 


(18c) 


Ineqns. (18), Ax and Ay are grid spacings along the x and y directions, respectively, while At is 
the time step. The six finite-difference equations are stepped in time, alternately updating the electric 


9 ^ 

and magnetic field components at each grid point. The value of o t j is about 10' ~10‘ and 


„ <7 At 

Ax = Ay = 10 m has been used here, so that: — « 1 . This condition simplifies eqn.(18a) to: 


i.j 


Er(iJ) = mj) , A t 1) 


«. j 


Ax 


Ay 


].(18d) 


This FDTD method is well suited for solving the cellular scattering problem since the 
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leapfrog method leads to a second-order accuracy in both time and space. Furthermore, the following 
Courant criterion as suggested for Yee’s algorithm was implemented to guarantee accuracy and 
stability [29] : 


v • At < • 

max — 


1 


1 


+- 


1 


lAx^ A/ Ar 


(19) 


where v max is the maximum signal phase velocity in the configuration being considered. If the value of 

one of Ax and At is known or assigned, the value of other one can then be obtained. 

Due to finiteness of the region used in any numerical simulation, care has to be taken to avoid 
artificial back-reflections from the simulation boundaries. Appropriate absorbing boundary 
conditions have to be implemented to ensure that unrealistic and unphysical wave reflections do not 
introduce spurious energy within the simulation region. A variety of boundary conditions have been 
reported [31] to prevent artificial reflections at the edges of the computational domain. These are 
either based on the solution of the wave equation [32,33], or the use of absorbing layers [30], The 
Berenger approach [30] of using a perfectly matched layer (PML) provides orders of magnitude 
improvement in performance relative to the other techniques [29], and so has been implemented here. 
Since the wave decays very rapidly in the PML material of onlyafew layers (—10), exponential time- 
stepping has been used instead of the usual Yee time-stepping. This changes equation ( 1 7) within the 
PML to: 


E n :\Uj) = e~ ff ^ ,e E n a ^ 

<J x (i) Ax 


[-H "(/ + 1/2, j) + //"(/- 1 / 2, y )] , (20a) 


E~'{Uj) = e-’ U) *"E m v - ' } ‘ e —— - [~H n x (/, y + 1 / 2) + H n x (i,j -l / 2)] , (20b) 

0 y (j )Ay 


H?'(i + l/2,j) = e 


-o\(j+V 2)A, >M H n ( /+ 1/2J) + 


a (j+l/2)At/ /i 


) 


cr y (y + l/2)Ay 


x {£"(/, j + 1) + E-JiJ + 1 )-[E'„(iJ) + El(i,j)]} 


(20c) 


. /I _ ~-<r'A‘+U *)*' ft \ 

H n+1 (f + 1 / 2, y) = e~° AM l2) *" u H n v (i + \/ 2, j)+ (l f } - 

<J X 0 + 1 / 2)Ax 

x {E n zx {i + l,y) + £;,(/ + l,y) ■ -[£-(/, j) + E n 2I {iJ)]} 


(20d) 
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In the above, a and cr* are the electric and magnetic conductivities, respectively, and both are 
functions of the grid position index i and j. Component E z is divided into Ezx and E zy in the PML for 
TM case. In the computation region, the traditional discretizations for the curl functions 18(a-c) are 
used. If the impedance of the PML medium equals that of the material used in the computation 
region, then no reflection occur at the interface of these two materials. An essential conditional 
requirement is: 

— = ^r ( 21 ) 

£cr Mcr 

The index "cr" denotes the computation region, and £ 0 and // 0 are used if the computation domain is 
air or a vacuum. 

One of the advantages of this time domain method is that both narrow- and broad-band source 
excitations can easily be obtained. For broad-band pulse signals, a discrete Fourier transform (DFT) 
of the transient behavior yields the complete frequency response. The discrete Fourier transform, 
taken at every time-step of the simulation, is expressed by eqn. (22) below: 

G(kLf) = A/5~’g(rtAf)exp( — 12^2. 1) } with k = 0,1,2, ...Nf , (22) 

to N, 

where G(kAf ) is the complex frequency domain data, g(nAt) are the time domain E and H fields, N t 
is the length of the DFT, A f is the frequency step, and Nf represents the number of frequencies 
N f =1/(A/A/). 

Due to the computational limitations, a relatively small simulation domain is chosen for 
obtaining the fields from the FDTD method. This region closely surrounds the cell, and so, only the 
near-field distributions are directly obtained. Far field information on the intensity and phase can be 
obtained from the near-field E and H values based on the free-space Green’s function and an 
integration procedure over a surface surrounding the cell [34,35]. In this procedure, the radiated E 
and H fields (in the frequency domain) tangential to a virtual surface completely surrounding the cell 
are converted into equivalent electric and magnetic surface current densities, J eq and Me q . These are 
then weighted by the Green’s function and finally integrated to provide the desired far-field response. 
The current densities J eq and M eq can be obtained as: 

J eq {r',co) = nxH{r'M (23a) 
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(23b) 


M eq {r',(d) = -hxE{r\(o) 

where r' is point on the virtual surface and h is a unit vector normal to the surface. H and E are 
obtained by averaging the nearest components around the grid point on the virtual surface. For 
example, in the 2D case, the far- field is produced by equivalent electric and magnetic currents that 
are tangential to the non-physical virtual contour as shown in Fig. 3. The figure shows the overall 
geometric schematic that can be used for the simulations. The cell is shown in the center, with the 
electromagnetic wave incident from the left. The PML surrounds the simulation region completely 
on the periphery. Two contours, one for the near-to-far field transformation, and the other for 
evaluations of the total and scattering fields have also been shown in Fig. 3. On the left edge of the 

virtual contour for the near-to-far field transformation, E(i + 1/2, j) and H(i, j ) are obtained by: 

E(i, j + 1 / 2) = [E(i, j ) + E(i, j + 1)] / 2 (24a) 

H(i + 1 / 2, j) = [H(i, j) + H(i + 1, j)] / 2 (24b) 

The scattered electrical field in the far field region is given by [36] : 



Fig. 3 Simple schematic of the geometry used for the FDTD scattering calculations for a single cell. 
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E s z (r\(0) = -jcoz'A 1 ( r\a>) + j(Mjr.F{r' ,oS) 

(25a) 

A z {r\(0)=^j[j s {r',(0)H^[k \r-r'\]dl' 

(25b) 

F z {r\(0) = ^[M s {r\(0)rH^[k \r-r'\]dl' 

(25c) 


Here, r is the point at far field. As | kr |— » °° , the Hankel function can be written as: 

H«\kr)~A vl e- jh . (25d) 

mr 

The phase function describing a scattering pattern, F((f>) is thus: 

F(<t>) ~ j== ^[-Tjz'-J eq (r')-zxM eq -r]e Jirr 'Ax . (25e) 

Here r\ = — , k is the wave number, N is the number of segment on the virtual contour, z is the 

V *0 

unit vector on z direction, r is the unit vector of the far field point and <p is the scattering angle. 

V. Validation of the FDTD Numerical Implementation 

Validation for the numerical implementation of the above FDTD mathematical model was first 
carried out by comparing the predicted scattering patterns for homogeneous spheres based on the 
FDTD technique, against the results of Mie theory [37]. The gird size for the FDTD calculations was 
X / 40 . The source pulse used in the simulations was a Gaussian pulse, propagating along the x 
direction defined as: 

E‘(s, j) = exp(- (t Jt - 3) 2 ) sin(2 #(/ - 3 r)) . (26) 

in eqn. ( 26 ), E K '(sj) is the line source, /is the central frequency of interest, and r is the characteristic time 
of the Gaussian pulse determined by the frequency bandwidth of interest. 

Scattering fields were obtained at every grid point, as the Gaussian pulse propagated through the 
entire computation region. The number of time steps for the DFT computation [the n in eqn.(6)] was 
chosen to be sufficiently large to ensure the pulse magnitude approached zero. The phase function 
F(0) was computed after the FDTD was completed, in 5-degree increments. Results of the 
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scattering patterns obtained from FDTD simulations for scatterers with two different indexes of 
refraction (m=l .02 and m=l .37) are shown in Figs. 4a and 4b. Results of the Mie theory are also 
shown for comparison. The radius of the cell (assumed to be a homogeneous sphere placed in 
vacuum) was taken to be 2 pm, while the chosen optical wavelength was also 2 pm. The plots of Fig. 
4 clearly show that the refraction index strongly affects the scattering pattern. The number of 
fluctuations increases with the refractive index. More importantly, results obtained from the FDTD 
technique match the predictions of Mie theory [37] fairly well. The agreement at small scattering 
angle (0~50°) is particularly good, with a slight deviation at higher angles. This may be due to 
reflection at the boundaries and dispersion of scattered light within the FDTD grid. Fig. 5 shows that 
the comparison between FDTD and Mie theory for a 0.535 pm circular object of dielectric constant 
1 .37, with extracellular fluid of relative permittivity 1 .35 surrounding it. The cellular values were 
chosen in keeping with published data [38,39], A 0.6328 pm wavelength was used for these 
calculations, since it corresponds to that of a He-Ne laser, an optical source that is commonly used. 



Fig. 4. Validation of the vacuum 2D TM FDTD code. The environment is vacuum. The index 
of refraction of cell is: (a) 1 .02, and (b)l .37. X =2 pm and a = 2 pm. 
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Fig.5. Validation of the 2D TM FDTD code. Index of refraction of cell is 1 .37 and that of the 
enviroment fluid is 1.35. He-Ne laser wavelength ( X =0.6328|im) is used, a =0.525(im. 

These results once again reveal close agreement between the FDTD results and the Mie theory, 
thereby validating the electromagnetic model and its numerical implementation. Furthermore, Fig. 5 
shows that the code is not only valid for simulations in vacuum, but also capable of treating spatially 
non-uniform relative permittivities. 

VI. The Alternating Direction Implicit (ADI) FDTD Method 

Though the FDTD method is widely used for simulating electromagnetic problems, it has one 
significant drawback. It is based on the Courant-Friedrich-Levy (CFL) condition [29] that places a 
constraint on the time step for a given grid spacing. Under the traditional FDTD, smaller grid 
spacing (as would be required for complicated geometries with inflexion points, or sharp comers, or 
sizes comparable to the wavelengths) translates into a smaller time step. This increases the 
computational time and burden significantly. Also, analyses for high-frequency signals would be 
slow and time-consuming. 

A method for overcoming this difficulty has recently been proposed [40] for two-dimensional 
problems that introduces an alternating direction implicit FDTD method (the ADI-FDTD approach). 
The method is based on the alternating direction implicit method of mathematics [41], and has been 
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applied to Yee’s staggered cell [28] to solve Maxwell’s equations. The ADI-FDTD scheme has 
recently been generalized to cover three-dimensional (3D) situations [42, 43]. 

Since the ADI-TDFD scheme is basically similar to the FDTD method already discussed, the 
details are not given here. However, the inherent equations that change due to the use of ADI, are 
presented and discussed below. The 3D space is discretized using the standard Yee-gird. Berenger’s 
perfectly matched layer (PML) is used around the 3D computation region as the absorption boundary 
condition. It requires field splitting. Then the 6 conventional Maxwell equations are changed into 12. 
ADI demands the one FDTD time step into two steps. The equations used in the first 1/2 At are: 


E -E 

W i+1/ 2 J, it xy i+l 1 2, j, k 


»+W2 «+w z H s 

“ ,+l/2,y+l/2,/t zy ,+i/2,y+l/2.* ' “ t+Ul,j-U2Ji V i+) 1 2.J-U 2,k * 


At/2 


Hz) 


e r /2 -e i 

n 

■»li+l/2.y.* «l 

1+3/2,;,* 

At 12 


>7b) 


|«+1 / 2 | 

/i 

£ -£J 


H/.y+i/2,* >*l 

(,;+]/ 2,* 

At/2 


lie) 


in+1/2 I 

1 

1 

t*1 


^liJ+1/2,* 

1 i,;+l / 2,* 

At/ 2 


lid) 


e r +1/2 -£ | 

» 


(,y,*+ 1/2 

At 12 


lie) 


in+1/2 I 

n 

1 


*y \ i,j t k+i / 2 

i,y,*+l/2 

At/2 


+ (j E 

y ** i+u 2,j,k 


+ <7 E 

2 » i+l / 2 ,j,k 


~ y* i,y+i / 2,Jt 


+ ^^.vx 


iJ+\ i 2,Jt 


+ G x i,j,k +\ / 2 


+ cr E 

y ijmu2 


Ay 


H + H — (H + H ) 

y* i+I/2,;,it+l/2 y* i+l / 2,;,*+l / 2 ^ - v * i+l/2,;,*-l/2 ^ t+1 1 2,;,*-l / 2 ' 

A Z 


n+1/2 , n+ l/2 _ " +1 ' 2 rj |" +,/2 x 

,,;+l/2,* + l/2 ”li,y+l/2,t+W2 K *> ij + l/2,k-\/2 ** U;+l /2,*-l / 2 7 

Az 

H \ n . jj | n —(H I" +H I” ) 

77 « I i+l / 2,7+1/ 2,k ^ 32 I i+l / 2,;+l / 2,k ' ^ li-1 / 2,;+l / 2,it « li-l/2.;+l/2.A * 

Ax 


£ fi-ri / 4 wtiu nTi/ 

// + // — (f-[ + // ) 

** i+l/2,/,Jt+l/2 yi i,;+l/2,*+l/2 ^ * ;,;+!/ 2,*-!/ 2 yz i,;+l/ 2>-l tl } 

Ax 


H, 


n i iy " _ / ir _i_ /y n \ 

I,;+l/2,A+l/2 ^ /,>+)/ 2,*+!/ 2 V ZJC 1,7-1/ 2.*+!/ 2 V i t >-I / 2,*+l / 2 7 

Ay 


(27f) 


M 


«+ 1/2 n 

^ i,y+I/2.*+l/2 ^ i,y+l/2.A+1/2 


A//2 


* A 

+ <j y H Xy i j +U2 k + )l2 


r I” i 17 n I " i If " ^ 

I <,/+!,*+!/ 2 ^ j,y+l,/t+l / 2 ^ “li,/,*+W2 ^ ijMVl' 


Ay 


(27g) 
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rr n+l/2 tt n 

& i,y+l/2,*+l/2 « r,y+l/ 2,*+l / 2 


At 12 

(27h) 


| n+l/2 

1 n 

H . 

-7/1 

1 *+1/ 2, y 4+1 / 2 • li+l/2,y\Jt+l/2 

A* 

At/2 

(27i) 


| n+l/2 

[n 

// 

-//J 

^li+l/2,y.*+l/2 ^li+l/2,y,*-+l/2 

A* 

At/2 

(27j) 


, _ i n+1 / 2 

„ ^ in 

// 

-H 

» 1 ,+i / 2,y+l / 2,* ** 1 i+l / 2,y + l / 2,* 

A* 

At/2 

(27k) 


1 n+1 / 2 

[ n 

// 

-//. 

™ 1 i+l / 2,y+l / 2,* zy \ i+l / 2,y+l / 2,k 


At/2 


+ < 7 * H 


xz li,y+l / 2,k+\ / 2 




i+l/2,y\*+l/2 


+ <J X H VX 

x y x i+l / 2,y',*+l / 2 


n+l/2 „ n+l/2 , „ n + l/2 _ n+l/2 

£* -j- ^ ^ ^ \ 

>* i,y+l / 2,£+l / 2 >* t,y+l / 2,*+l / 2 V v * i,y+l / 2,*-l / 2 >* i,y+l / 2,*-!/ 2 7 

Az 


.'+I/2.M+I ^ x >’ MH.j.k +£ ”1 i> 1/2,;,*^ 

Az 


n+l/2 n + l/2 _ n+i/2 " +1/2 v 

^ i+i,;,*+i/2 M,jM\u ” I.j*+\n *y IJMU2/ 

Ax 


+ a r H 


• TT I" 


I i + l/ 2,y+l/ 2,k 


E j? (E ^ \ 

y* i+l,y'+l / 2,k yz t + l,y+l / 2,k ' yX i,j+\/2,k >* li,y+l / 2,;T 

Ax 


+ (7 H 

y *y i+l/2,y+l/2,it 


n+1 / 2 

^ i+l/2,y+U 


, F " +1/2 _*£ " +l/J ,£ n+,/2 \ 

Tc, xz i+l / 2,j+l,k ' » 1+1/2, ;.* O' f+i/2,;,* ' 


Ay 


(27m) 


The equations used in the second 1/2 At are: 


n+1 n+l/2 

/T — F 

*>’ i+l / 2,y,Jt ^ i’+l / 2,y\Ar „ |«+^2 

At/2 y 


n+ 1/2 n+l/2 

— ( H n ' * 4- // 1 

v 2 * zx , + ) ( 2,y- 1 ' ’ * n : i i l 2 


f~f ” + *^ i J-T 

i+l / 2,2+1/ 2.* O' ,+1/ 2,2+1/ 2,/t v " a li+l/2.2-l/2.2 ' " O’ I/+I /2.2-I/ 2.* 


(28a) 


| n+1 


xz li+l/2,y,* 


_ n+l/2 

J.k ^ xz i+\nj t k 


At/2 


n+l/2 

+ <X z^xi /+!/ 2,2.* 


Ay 


n+1 n+1 n+1 „„ n + 1 

^ _l_ {hi. 4 H ) 

y* i+l / 2,j,k+V2 y* 1+1/2,;,*+) / 2 ^ yx i+l/2,y,*-l/2 yz i+l / 2,y,*-l / 2 ' 

Az 


(28b) 


n+1 n+l/2 

E -E 

y 2 i,/+l/2,Jt y* rJ+l/2,* | l »+l / 2 

£ - + <7 E 

At/ 2 2 yr| ^ +1/2 ’* 


n+l/2 n+l/2 _ f Z + // * +1 ' 2 A 

*y ij+ 1/2, *+1/2 ** «,;>!/ 2,^+1/ 2 V *> i,y+l/2,*-l/2 ” /,y+l/ 2,*-!/ 2 ^ 


n+l/2 


Az 
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(28c) 


in+1/2 


7 

' yx tj+\/ 2 ,k yx ij+\n,k 


At/2 


n+ 1/2 

+<T * £ y « (,y+i/ 2 ,* 


H 


ln+1 


+H„ 


— ( hf i J-f 'i 

«li+l/2,y+l/2,* ' “«li+l/2,y+l/2,* ,-l/2,y+l/2 ( * « M/2,y+l/2X ; 

Ax 


(28d) 


| M+1 /i n+1/2 n+1/2 n+1/2 n+1/2 

n I ri I 2 ^ _j_ ^ ^ ^ j 

^ or I ,• j t *+l / 2 ^“liJX+1/2 _ jy | n+1 / 2 >* 1+1/ 2,7,*+]/ 2 >* (\y+l/2,*+l/2 ^ >* i,y+l/2,*-l/2 1,7+1 1 2,*-l 11 ' 

£ — + <7 * E *\,JMU2 ~ 


At / 2 


(28e) 


I n+1 / 2 


rt + 1 ft ^ 

r _ 17 

^ »,y,*+l/2 ^ i,y,*+l/2 ~ |«+l/2 

£ Z771 + a ’ E X,> ..<> 

(28f) 

in +1 in+ 1/2 

■ l H/.;+l/ 2 .*+l /2 J y|i.j>l/ 2 ,t-H /2 ^ \n+Vl 

" At/ 2 y J °'l>,;+i/ 2 ,/t+i /2 


n+ l " +l n+1 " +l x 

“ /,;+l/2,*+l/2 V (,;+l/2.*+l/2 ' r * .,y— 1 1 2,*+l / 2 ^ ,,;-|< 2,*+l /2 _ 


Av 


^ «+l r " +l " +1 ,£■ ' I+1 \ 

** i\y+l,*+l/2 *>’ /,y+l,* + l/2 ^ zx iJMVl *>’ i.y'X+l / 2 


Ay 


( 28 g) 


__ n+1 jj n+1/ 2 

i,y+l/2,*+l/2 ~~ « i,y+l/2,*+l/2 , I n+1/ 2 

V » . / “ «|,,y + l/2jl+l/2 


At! 2 


(28h) 


n+1 n+1/ 2 

// — // 

^ t+l/2 t ;,*+l/2 i+W 2 ,j ,*+1/2 ^ q* J-[ |” +1/2 

^ A//2 1 i+l/ 2,;,*+!/ 2 


n+1/ 2 n+1/ 2 n+1/ 2 n+1/ 2 

y* i tj / + l / 2,*+l / 2 >* i,y+l/2,*+l/2 ^ yx i,y+l / 2,*-l / 2 - vz iyilWT 

A Z 


7 " +1 r «+l 0 + 1 n+1 

/+ 1 / 2 .M +1 +tja ,+ 1 / 2 *M+ l /+ 1 / 2 .M _ * z 

Az 


(28i) 


n+1 n+1/2 

^ — // 

>* i+l/ 2,/,Jt+l / 2 y x i+l/2,y,*+I/2 , * JJ 

U + 

• A X / ^ * >* 


£ .r + . 1/2 . .... +£’. 1,1+1 2 


n+1/2 n+1/2 

-( £ -U.,. +£ «l, 


41/2 


n+1/2 _ z*Ii+1,7,*+1/2 ^Ii+1,»+1/2 

i+ 1/ 2,y,Jfc+l/ 2 A* 


(28j) 


rr « +1 _zy 0+1/2 

11 ** i+l/ 2,7+1 /2X « i+l/ 2,7 + 1/ 2.* , rr I n+1/ 2 

M : — + 0 ^ 


A// 2 


X »l(+l/ 2 , 7+1 / 2 ,* 


(28k) 


n+1 n+1/2 

// — // 

^y t+i/ 2,7+1/ 2,* z y t+i/ 2,7+1 11M | in+1/2 

^ A// 2 ^ ^ •»+! / 2,7+1 / 2 t k 


n+1 n+1 n+1 n+1 

^y* i+1 ,7+1/2> + - vz 1 + 1,7 + i /2,t ^ y j 1,7+1/ 2,* + > T »,7+1/2X^ 

Ax 


I n+1 2 I n + 1 / 2 

+£ .|, 

Ay 


<r+l/ 2 , r< n+1/2 

/+i/2,/+i,t ” i+i/2,;+i,* v " ”11+1/2,;,* ’ ~^l(+i/2,;,* 


(28m) 
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Bring (27k) and (27m) into (27a), we get: 

- M + i/2 jjk BB y\ M <2,H,k + (1 + By \ {BB ^ 


i+ 1 / 2J,k 


BBy | 


n+1/2 

i+l/2,7-1,*^^ |+|/ 2 , y, it 


-^Ll/2,y.* 55 -H W/ 2,y,* £ 

Here H i +,/2J + l/2 > *=^ 


rt + 1/2 


= Mi 


M/2J+Uk 
n+1/2 | n+1/2 

i+l / 2,y +1 / 2,* ** I i+l / 2J~\ ! 2,k 


zj 

i+l/2J+l/2,k ** 


*'+l -2J t k *y 1 1+ 1/2 ,y 

+ Kai Lv2J + V2,k~ Ka Lu2. 

+ B H + u 


(27n) 


i+l/2./+l/2Jt 


li+l/2,y-l/2,A 7 


2./+l/2.Jt v “I/./+U 


(27n) is the implicit update expression of Exy, a tri-diagonal matrix needs to be solved. In the 
same way, we can get the update expression of Eyz and Ezx in the first half step. 

Bring (28i) and (28j) into (28a), we get: 




li+l/2,;,*-l 


+ 0 + *U« + BBz I, 


- Bz BBz E 

*1/2 J,* *1/2 JJc x 




in+1/2 


n+1/2 

i+l/2,y,*-l 

n+1/2 I I" 

M/2JMl = AZ ' MI2 'J' k 

n+1/2 


n+1/2 
i+l/ 2j\k 


(28n) 


1+1/2,/,* v y x |i+l/2,/,*+l/2 


— H 


yx \i*V2,j,k-\l2 


+ Kai\ 


i+l/ 2 ,y,*+l /2 


-H 


i+l ( 2J t k-\ ( 2 J 


Here&7i|. 


i+l /2,_/,^+l/ 2 


= AAzl 7/ " + 55z 

li+l/2,y,*+l/2 /+l/2,y,*+l/2 


,(*, 


j n+1/2 


i+l/ 2 ,y,*+l/ 2 v ^li+l/2J,*+l 




n+1/2 
i+l / 2,y\* 


) 


In the same way, we can get the update expression of Eyx and Ezy in the second half step. 
The coefficient used in (27n) and (28n) are: 


M,jjt = 


1 + 


°jk&_ 

jf £ 

4 fy* 


,q = x,y,z , 
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H„ = 


At 

2e iik Aq 


1 + - 


'V k 

VjpA<_ 


,q = x,y,z , 


1 — 




=' 


1 + 


, ,q = x,y,z , and 
cr„,At 


ij* 

Afi 


ijk 


At 


BBq \ijk = 2M '^l ,q = x,y, 
C ijk At 


1 + - 


4// 


ijk 


VII. Simple Example Implementation and Results 

Numerical codes were developed and implemented for electromagnetic calculations of the 
scattered fields and power, based on both the MM and hybrid schemes. The intent was to 
demonstrate that the simulation capability is successfully in place at Old Dominion University. These 
codes could be tested further for validation and refinement, and applied to real scattering problems 
with complex geometries. Here, the numerical codes were used only for a simple geometry and the 
calculation results presented in detail in the next section. 

VII.A MM - Related Calculations 

A simple flat, two-dimensional (2D) rectangular geometry was considered lying in the z = 0 
plane. A frequency of 30 GHz was chosen. In accordance with the triangular patch basis functions 
for surfaces given in section III. A, this 2D plate was divided into 16 triangles with uniform grid 
spacings Ax = Ay (s A) , which was taken to be 1/10 of the operating wavelength. The extent along 
the x-axis was 4 A, and 2A along the y-axis. In general, rectangular shapes with M and N segments 
along the x- and y-axes lead to 2MN triangles, with a total number of basis functions Nb as is (which 
requires an accounting of all mutually adjacent triangles) given by : Nbasis = 3MN - M - N. Hence, 
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this geometry required 1 8 basis functions. All triangles had identical areas. 


The numerical implementation has been carried out completely, with the corresponding 
source code written in FORTRAN. Results for the impedance matrix [Z] and the current were 


obtained, and are given below. The [Z] elements, which were complex, were as: 


1 1 (3. 4616823 E- 6,1. 284088 6E-3) 

1 3 (3 .0042051E-6, 3 . 5189796E-6) 

1 5 (3. 842653 E- 6, -4. 765277 5E-4) 

1 7 (3. 645482 E- 6, 5.2011955E-6) 

1 9 (6.489668E-6,3.6305928E-6) 

1 11 (4.454643E-6,3.6243454E-7) 

1 13 (7 .99366E-6, 4.103712E-5) 

1 15 (5.052641E-6, 6.984539E-7) 

1 17 (5. 4896554 E- 7, -3. 2988587 E- 6) 

2 1 (3.4260583E-6,2.027901E-5) 

2 3 (3.4333897E-6,2.0456453E-5) 

2 5 ( 3 . 994305E-6 , -4.872825E-4) 

2 7 (4. 065163 E- 6, 2. 2113232 E- 5) 

2 9 (7 . 05314 05E-6 , -2.9536808E-5) 

2 11 (6.072316E-6, 4 .264207E-6) 

2 13 (8. 248924 E- 6, 5. 155328 E- 4) 

2 15 (6. 9864927 E- 6, 7.718715E-6) 

2 17 (8.0898826E-7, -1. 2277626E-5) 

3 1 (3.0029814E-6,3.5050061E-6) 

3 3 (3.4605432E-6,1.2235348E-3) 

3 5 (3 .7368145E-6, 2 .0052641E-5) 

3 7 (4. 0803574 E- 6, -4.995696E-4) 

3 9 ( 6 . 8089975E-6 ,6.1009304E-6) 

3 11 (7.126019E-6,7 . 5248003 E- 6) 

3 13 (7. 513592 6E-6, 3.8112284E-6) 

3 15 (8.228179E-6,4 .514658E-5) 

3 17 (1. 068182 E- 6, 2 . 3982657 E- 5) 

4 1 (2. 2702647 E- 6, 5. 3112523 E- 7) 

4 3 (3 . 0962701E-6, 1.5218275E-5) 

4 5 (3.179949E-6,4.42893E-6) 

4 7 (3. 75 57 738 E- 6, -5.070338E-4) 

4 9 ( 5 . 7958486E-6 , 3 . 9000697 E- 6) 

4 11 (7 . 384672 E- 6, -2 . 3749117E-5) 

4 13 (5.9176313E-6, 1.258488E-6) 

4 15 (8.487027E-6,4.922489E-4) 

4 17 (1. 2615717 E- 6, 1.427646E-5) 

5 1 ( 1 . 3282961E-6 , -5.065996E-4) 

5 3 (1.5254552E-6,1.7916516E-5) 

5 5 ( -3 . 390591E-6 , 5.512347E-4) 

5 7 (-1. 4370028 E- 6,7. 5093394E-6) 

5 9 (-1.4428904E-6, 2 . 5558926E-5) 

5 11 (1.0101412E-6, 6.0330712E-6) 

5 13 (-1.857216E-7,-3.6293978E-4) 


1 2 ( 3 . 2690945E-6 , 1. 7910324 E- 5) 

1 4 (2.624891E-6,8.7982584E-7) 

1 6 (3. 7359914 E- 6, 2. 01 80203 E- 5) 

1 8 ( 6 . 8873705E-6 , -3. 243064 E- 5) 

1 10 ( 5 . 65048E-6 ,3. 3914564 E- 6) 

1 12 (8.129872E-6,5.271748E-4) 

1 14 (6 . 8798726E-6, 7.3656983E-6) 

1 16 ( 8 . 0954805E-7 ,-1.2388141E-5) 

1 18 ( 6 . 7900316E-7 ,-2.863137E-7) 

2 2 (3.4611057E-6,1.2538116E-3) 

2 4 ( 3 . 2 9 587 2 E- 6 , 4. 305291 E- 6) 

2 6 (3.9615124E-6,-4.8804885E-4) 

2 8 (6.809017E-6,6.1009614E-6) 

2 10 ( 6 . 8078406E-6 , 5 . 577691E-6 ) 

2 12 (7. 40168 E-6, 3. 0495348 E- 6) 

2 14 (8. 110917 E-6, 4. 3091848 E- 5) 

2 16 (1.2968917E-6,2.6661093E-5) 

2 18 ( 7 . 586889E-7 ,-2.643961lE-6) 

3 2 (3. 2 611687 E- 6,1. 7748624 E- 5) 

3 4 ( 3 . 5976945E-6 , 2.3002576E-5) 

3 6 (3.875039E-6,-4.971584E-4) 

3 8 ( 5 . 9 3 647 IE- 6 , 4. 206562 E-6) 

3 10 (7 . 2189095E-6, -2.6642956E-5) 

3 12 (5.821322E-6,1.0841968E-6) 

3 14 (8 .367975E-6, 5.038911E-4) 

3 16 (1. 6828 57 4E- 6,1. 5951882 E- 5) 

3 18 (8 .084246E-7, -1 . 2167126E-5 ) 

4 2 (2.7109927E-6,2.7208943E-6) 

4 4 (3.4599738E-6,1.193258E-3) 

4 6 ( 3 . 4 06967 IE- 6 , 1 . 812755E-5 ) 

4 8 (4 . 4610433E-6 , 4 . 2454996E-7 ) 

4 10 (6. 8089884 E-6, 6. 1009163 E-6) 

4 12 (3 . 7453246E-6, -1.5543089E-6) 

4 14 (7. 62 54 937 E-6, 4. 572916E-6) 

4 16 (1.8918652E-6,4. 909446E-6) 

4 18 (8 . 3947156E-7, 2.1304206E-5) 

5 2 (1. 2004091E-6, -5. 121067 E-4) 

5 4 (2. 11892 57 E-6, 5. 4923284 E-6) 

5 6 (-2.9263424E-6,2.1411299E-5) 

5 8 (-1 . 8144461E-6, -1. 8135617 E- 5) 

5 10 (-3.9065594E-7,1.3647906E-5) 

5 12 (-3.9784663E-7,-4.2573123E-5) 
5 14 (-2.5795543E-9,6.418809E-6) 
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5 15 ( 

5 17 ( 

6 1 (1 
6 3 (1 
6 5 (- 
6 7 ( - 
6 9 ( - 
6 11 ( 
6 13 ( 
6 15 ( 

6 17 ( 

7 1 (1 
7 3 (1 
7 5 (- 
7 7 (- 
7 9 (- 
7 11 ( 
7 13 ( 
7 15 ( 

7 17 ( 

8 1 ( - 
8 3 (5 
8 5 (- 
8 7 (- 
8 9 ( - 
8 11 ( 
8 13 ( 
8 15 ( 

8 17 ( 

9 1 ( - 
9 3 (2 
9 5 (- 
9 7 ( - 
9 9 ( - 
9 11 ( 
9 13 ( 
9 15 ( 

9 17 ( 

10 1 ( 
10 3 ( 
10 5 ( 
10 7 ( 
10 9 ( 
10 11 
10 13 
10 15 

10 17 

11 1 ( 
11 3 ( 
11 5 ( 
11 7 ( 
11 9 ( 


1.2666419E-7, 8. 195884 E- 7) 

-9. 182856 5E -7,1. 8520703 E- 6) 

.. 5312934E-6, 1 . 7810573E-5) 

.. 3194917 E- 6,-5. 397907 E- 4) 

2 .9398602E-6, 2. 1336144 E- 5) 

2. 6039947 E- 6, 2. 6274956E-5) 

1 . 698 9104E-6 , -1 . 696899E-5 ) 

1. 3675981 E- 7,1. 57 00189E-5) 

-2 . 3212505E-7 , -3 . 8739046E-5 ) 
1 . 54084 55E-7 , 7 .4791037E-6) 
-1. 4137604 E- 6, -3.378447E-5) 

. 6527874 E- 6, 4 . 688458 E- 6) 

. 08797 46E-6 , -5. 368387 E- 4) 

2. 0019865 E- 6, 5. 8753985 E- 6) 

3. 39273 IE -6,4. 90687 55E- 4 ) 

1 . 62 55885E-6 , -4 . 2 83884E-6 ) 
-7. 550989E-7, 3.4162192E-5) 

-3 . 00342 54E-7 , - 5 . 29713 72E-6 ) 
1. 5023443E-7 , -3 . 9793216E-4 ) 
-1. 5414896E-6, -3 . 4250051E-6 ) 
2.4519601E-7, -4. 0247177 E- 5) 
.7997033E-7, 3 .2883967E-6) 

4 . 2161655E-6 , -1 . 7798244E-5 ) 
2. 171954 E- 6, 5. 6421765 E- 7) 
3.8866637E-6, 9. 593817 E- 6) 

-9. 785879 E- 7, 2 .4496607E-6) 

-2 . 95137 E- 6, -3 . 568135E-5) 

-2 . 1892651E-6 , -1.5354407E-6) 
-6. 315448 E- 7,1. 8196948 E- 5) 

1 . 47 9827 5E-8 ,1. 2888796 E- 7) 

. 97 62952E-7,1. 854493 6E-6) 

3. 972510 5E -6, 1.6900563E-5) 
3.3975948E-6, -3. 9080623 E- 6) 
4.3537E-6,1.1814693E-3) 

-2.3 3 61235E-6 , 3 .0384974E-6) 
-2. 997729E-6, - 5 . 058632E-4 ) 
-2. 8159452 E- 6, -6.4250875E-6) 
-1.2 655715E-6, - 5 . 1834 876E-4 ) 
3. 0580008E-7,2.8626718E-6) 

8. 21002 E- 8, -3 .449577E-5) 
-3.123678E-6, 1.162 57 51E-5) 
-4. 214612 E- 6, -1. 7551507 E- 5) 
-4. 206985 E- 6, 4. 17734 E- 6) 

(-3 . 5572693E-6,1. 4693527 E- 5) 
(-2 .699236E-6, 5.02313E-7) 

(-3 .186187E-6, -3 .9790753E-5) 
(-1. 6546347 E- 6, -5.261443E-4) 
6.2651E-7,1.7669165E-6) 

-1. 8407604 E- 8,1. 0251813 E- 7) 
-1. 8757693 E- 6,4. 3743116 E- 6) 
-4 . 4270946E-6 , 1.156214E-5) 
-3.498722E-6, -1.2194391E-7) 


5 16 (-1.2482839E-6, -3.116474E-5) 

5 18 (2.233487E-7,5.4667475E-6) 

6 2 (1. 208135 5E-6, -5. 2171933 E- 4) 

6 4 (1.8564879E-6,2.0057238E-5) 

6 6 (-3.3916628E-6,5. 2096135E-4 ) 

6 8 (-1.5219584E-6, -4. 0074432 E- 6) 

6 10 (-1.0989961E-6/2. 9860565 E- 5) 

6 12 (-4.518921E-7,-5.8954824E-6) 

6 14 (-1.774208lE-8,-3.8043613E-4) 

6 16 (-1.0718904E-6,-5.079302E-7) 

6 18 (-7. 662204 E-7, 2. 6550505 E- 6) 

7 2 (1. 1994189 E- 6,1. 5679585E-5 ) 

7 4 ( 1 . 4385755E-6 , -5. 6747434 E- 4) 

7 6 (-3.2622435E-6,1.6494674E-5) 

7 8 (-9. 989093 E-7, -6.683132E-8) 

7 10 (-1. 5833828 E- 6, -1. 5802372 E- 5) 

7 12 (-4.809328E-7, -1. 7230377 E- 6) 

7 14 (-6.639891E-8, -3. 490495 E- 5) 

7 16 (-5.961662E-7,4. 2460965E-6 ) 

7 18 (-1.5792373E-6,-3.640419E-5) 

8 2 (-1.8736351E-8, -8.539237E-8) 

8 4 (1.2973869E-6, 2 . 293018 E- 6) 

8 6 (-3. 608613 E- 6, -4.5680812E-6) 

8 8 (-4.3539766E-6,1.2117194E-3) 

8 10 (-2.6280004E-6,2.2516192E-6) 

8 12 (-2 .8783417E-6, -4. 96115 IE -4) 

8 14 (-2.7093383E-6,-6.0720803E-6) 

8 16 (-1.3838471E-6,-5.068088E-4) 

8 18 (6.16491E-7,5.833969E-6) 

9 2 (-8. 154793 E- 8, -3. 7371471 E- 5) 

9 4 (1.000247E-6,4.158571E-6) 

9 6 (-4. 215389 E- 6,-1. 7674862 E- 5) 

9 8 (-4.04241E-6,6.700779E-6) 

9 10 (-3. 721967 E- 6,1. 2143675E-5 ) 

9 12 (-2. 5874437 E- 6,1. 263982 5E-6 ) 

9 14 (-3. 0687801 E- 6, -3 .7736066E-5) 

9 16 (-1.5348432E-6,-5.162637E-4) 

9 18 (-3.0247724E-7,2.0127358E-5) 

10 2 (-1.6605337E-8, 1.1569311 E-7) 

10 4 (6.139951E-7, 3. 7944073 E- 6) 

10 6 (-4. 199800 5E-6,1.4231362E-5) 

10 8 (-3 .2069831E-6, 6. 613143 E-7) 

10 10 (-4.353419E-6,1.151219lE-3) 

10 12 (-2.0479958E-6, -5. 5485156 E-7) 
10 14 (-3 . 1171157 E- 6, -5. 1561126E-4 ) 
10 16 (-1.2910064E-6,1.4239511E-5) 

10 18 (-1.1472957E-6,-5.298888E-4) 

11 2 (1.6354984E-7,2.5524761E-6) 

11 4 (2.4574524E-7,-3.1620053E-5) 

11 6 ( - 3 . 543 8156E-6 , 9.9529424E-6) 

11 8 (-2. 036036 6E-6,1. 4151904 E- 6) 

11 10 (-4. 3715613 E- 6,1. 6539116 E- 6) 
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11 11 (-4. 353149E- 6,1. 1209687 E-3) 

11 13 (-2 .1443024E-6, -7 .2913536E-7) 
11 15 (-3.236504lE-6,-5.2535894E-4) 

11 17 (-1. 620924 6E-6,1.2313833E-5) 

12 1 (4.86008E-6,5.2826735E-4) 

12 3 (3 .4827353E-6, 1. 0371624E-6) 

12 5 (4.0270933E-6,-3.307024E-5) 

12 7 (2.268097E-6,-2.0473858E-6) 

12 9 (5.07908E-6, 5. 9779922 E- 6) 

12 11 (2.0555075E-6,-2.1821397E-6) 
12 13 (4.359218E-6,-1.3717479E-5) 

12 15 (2.1456116E-6,-2.2294877E-6) 

12 17 (2.502673lE-6,6.01409lE-6) 

13 1 (4.7143912E-6,3.7682817E-5) 

13 3 (4 .371323E-6, 9 .165892E-7) 

13 5 (4.211872E-6,-3.24895E-4) 

13 7 (3.3735E-6, -4.061778E-6) 

13 9 (5. 808194 E- 6, -4.4598E-4) 

13 11 (3 .7652316E-6, -7 . 6195806E-7) 
13 13 (4.6897002E-6,6.35869E-4) 

13 15 (3.505615E-6,-2.2868644E-6) 

13 17 (3.1629147E-6,4.0453294E-5) 

14 1 (4.045937E-6,6.0831753E-6) 

14 3 (4.7829576E-6, 4. 95624 E- 4) 

14 5 (3.7910394E-6, 8.388228E-6) 

14 7 (4. 084074 E- 6, -3. 038997 E- 5) 

14 9 (5.7168763E-6, -1. 932701E-5) 

14 11 (5.1499346E-6, 5. 800220 6E-6) 

14 13 (4. 3222653 E- 6,-1. 06866 5E-5) 

14 15 (4 .469397E-6, -1.23309E-5) 

14 17 ( 3 . 507475E-6 , 3. 679107 E- 4) 

15 1 (2.9906414E-6.1.0069923E-6) 

15 3 (4.6364134E-6,3.566128E-5) 

15 5 (2. 8951822 E- 6, 6.0595147E-7) 

15 7 (4 . 2640926E-6 ,-3.3341735E-4) 

15 9 (4.823016E-6,-1.7244968E-6) 

15 11 (5.886068E-6, -4 . 133595E-4 ) 

15 13 (3.3484016E-6,-2.0766692E-6) 
15 15 (4. 803733 E- 6,6. 267012 E- 4) 

15 17 (3. 419923 5E-6,1.1537359E-6) 

16 1 (3. 1832396E- 6,-1. 2917593 E- 5) 

16 3 (3.746049E-6,1. 5019422 E- 5) 

16 5 (2.4329631E-6,-2.4869798E-5) 
167 (2.7302235E-6, 5. 62648 E- 6) 

16 9 (6. 5537283 E- 6, -4.780589E-4) 
1611 (5.401115E-6, 5.2356112E-6) 

16 13 (7. 816293 E- 6,3. 265896 5E-4) 

16 15 (5.09427E-6,-1.9881284E-7) 

16 17 (4.943645E-7,2.4520778E-5) 

17 1 (2.9264238E-6,-1.7075954E-6) 

17 3 (3.7765526E-6, 3. 432097 E- 5) 

17 5 (2.52383E-6,3.9494016E-6) 


11 12 (-1 . 3598656E-6 , 3.383081E-7) 

11 14 (-2 . 8110275E-6, -2 . 593456E-7 ) 
11 16 (-7. 77504 8E- 7, 3. 555 57 72 E- 6) 

11 18 ( -1 . 7744276E-6 , -5. 3602457 E- 4) 

12 2 (4 .4065127E-6, 8.269442E-7) 

12 4 (2 . 2930721E-6, -6.2110854E-7) 

12 6 (3.3456917E-6,-4.278911E-6) 

12 8 ( 5 . 769255E-6 , -4.6229013E-4) 

12 10 ( 3 . 7358036E-6 , -7 . 795356E-7 ) 

12 12 ( 4 . 6326836E-6 , 6. 4045283 E- 4) 

12 14 (3.4564918E-6,-2.4191292E-6) 
12 16 ( 3 . 1911517E-6 , 4.1785286E-5) 

12 18 (1.6217407E-6.7.148078E-7) 

13 2 ( 4 . 8215156E-6 , 5.119455E-4) 

13 4 (3 .4536001E-6, 1. 0199728 E- 6) 

13 6 ( 4 . 0555864E-6, -3.173011E-5) 

13 8 (5.6774297E-6,-2.0338608E-5) 

13 10 ( 5 . 1145066E-6 , 5. 8891073 E- 6) 

13 12 ( 4 . 267819E-6 , -1. 1276143 E- 5) 

13 14 ( 4 . 414308E-6 , -1. 3024182 E- 5) 

13 16 (3.5342582E-6,3.63650lE-4) 

13 18 ( 2 . 4745818E-6 , 5.796431E-6) 

14 2 ( 4 . 675405E-6, 3.6672034E-5) 

14 4 (4 .3361346E-6, 1.0062317E-6) 

14 6 ( 4 . 237981E-6 , -3. 291563 IE-4) 

14 8 ( 4 . 7860794E-6 , -1 . 8788076E-6 ) 

14 10 ( 5 . 8471296E-6 , -4. 2966983 E- 4) 
14 12 (3.3006272E-6.-2.1873131E-6) 
14 14 ( 4 . 7467165E-6 , 6 . 3128536E-4 ) 

14 16 (3.442733E-6,1.0234307E-6) 

14 18 (3 .1346801E-6, 3. 9121303 E- 5) 

15 2 (4.0093536E-6, 5 . 9297226E-6) 

15 4 (4.7443963E-6,4.793020lE-4) 

15 6 ( 3 . 8135476E-6 , 8.257315E-6) 

15 8 ( 3 . 3111055E-6 , -2. 0225952 E- 6) 

15 10 (5.756332E-6,-1.8315402E-5) 

15 12 ( 1 . 9609688E-6 , -2.219029E-6) 

15 14 (4 .3767076E-6, -1.00971665E-5) 
15 16 (2 . 9495506E-6 , 1.1702908E-6) 

15 18 (3.4806917E-6,3.7217108E-4) 

16 2 (3. 435 67 E- 6, 3. 0039517 E- 5) 

16 4 ( 3 . 9389565E-6 , 5.3375906E-6) 

16 6 (2.3665334E-6,4.399628E-6) 

16 8 ( 6 . 785991E-6 , -4. 7987667 E- 4) 

16 10 ( 6 . 0897586E-6 , 2.0778683E-5) 

16 12 (7 .867391E-6, 5.731445E-5) 

16 14 (6.8187973E-6, 6. 302161 6E-8) 

16 16 (4.5478827E-7,5.852904E-4) 

16 18 ( 1 . 1977281E-6 , 7. 844572 E- 6) 

17 2 (3.2961843E-6,-1.1774497E-5) 

17 4 (4.271087E-6,1.7065055E-5) 

17 6 ( 2 . 2696715E-6 , -2.7472983E-5) 
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17 7 (2. 520416 5E-6, 5.2169916E-6) 

17 9 (6. 666465 6E-6, -4.949858E-4) 

17 11 (6.4213836E-6,2.2922446E-5) 
17 13 (7.701049E-6, 5. 3480434 E- 5) 

17 15 ( 6 . 66197 9 5E-6 ,-9.972264E-7) 

17 17 ( 4 . 5446495E-7 ,5. 5500923 E- 4) 

18 1 (2.3686616E-6,-2.2426502E-7) 
18 3 (3. 4091262 E- 6,-1. 0631396E- 5) 
18 5 (2.4565272E-6, 5.6154176E-6) 

18 7 (2.1063815E-6,-3.007615E-5) 

18 9 (6. 0832644 E- 6, 2. 0650237 E- 5) 

18 11 (6. 7933647 E- 6, -5. 3334643 E- 4) 
18 13 (6.64835E-6,8.959851E-6) 

18 15 (7 . 4807785E-6 ,3. 0108 17 4 E- 4) 
18 17 (5. 139923 E- 7,2. 4796653 E- 5) 

The current magnitudes were as: 


1 (-7.0488726E-4,-4.1359872E-3) 

3 (1.1044028E-3,9.498428E-3) 

5 (-1.6145814E-3,-7.883782E-3) 

7 (1.9769833E-3, 0.015651448) 

9 (-9.280766E-4, -2 . 6803847E-3) 
11 (6.84211lE-4,5.6747663E-3) 

13 (-6.516619E-4,-2.5269184E-3) 
15 (-9.80073E-4,-3.3341408E-3) 
17 (-1.5530751E-4,3.6133027E-3) 


17 8 (6 . 4145865E-6, 2. 2780714 E- 5) 
17 10 ( 6 . 6735424E-6, -5. 057028 E-4) 
17 12 (6.799824E-6,9.558169E-6) 

17 14 (7.648538E-6,3.1383574E-4) 
17 16 (8.3751774E-7,2.9638471E-5) 

17 18 ( 8 . 155057E-7 ,2. 9356284 E- 5) 

18 2 (2.820383E-6,-1.993386E-6) 

18 4 (4 . 117 4 4 IE- 6 ,3.8602403E-5) 

18 6 (2 . 055875E-6 ,1. 0388953 E- 6) 

18 8 (5.3973044E-6,5.2112973E-6) 
18 10 (6.546939E-6, -5.1009486E-4) 
18 12 ( 5 . 0371545E-6 , 1.2531445E-6) 
18 14 (7 . 5347065E-6, 4. 9646412 E- 5) 
18 16 ( 1 . 2029079E-6 , 7.8511675E-6) 
18 18 (4. 54 13707 E-7, 5. 2472803 E-4) 


2 ( -3 . 259895E-4 ,3.7869822E-4) 

4 (1 . 08 5154 IE- 3 ,7.1884454E-3) 

6 ( 4 . 3422586E-4 , 7. 928428 E- 3) 

8 (-6. 5302517 E-4, -2. 8714612 E- 3) 
10 (8.852718E-4,9.267416E-3) 

12 (-1.2921324E-5,1.349302E-3) 

14 (-3. 27375 54E-4, 5. 139351 E-4) 

16 (-1. 3548372 E- 3, -6. 1254217 E- 3) 
18 (2. 0707286E-3,0. 016081538) . 


Code for the hybrid MoM-PO technique was similarly developed. However, the hybrid 
MoM-PO scheme simplifies for a planar geometry for the following reason. In the hybrid scheme, 
the full, self-consistent solution is to be obtained through an iterative procedure. First, the currents in 
the PO region can be taken to be zero (an initial guess), while the usual MoM solution applied to the 
MoM zone. This yields the current distributions Ji over the MoM region by solving : L e [Ji(r)] + 
L e [J 2 (r)] = - E lnC | Un (r) , for r e MMZ , and J 2 (r) = 0. Next, this value of Ji is used to update 

the current distribution in the PO region according to: J 2 (r) = 2 n A x H lnc (r) + Lh[Ji(r)] + 


Lh[J 2 (r)] , for r e POSR , where n A is the unit surface normal, and the operators L e [J(r)] and 

Lh[J(r)] defined previously in eqn. (14). For planar structure, both J(r’) and V G(r,r’) are in the 
plane. Hence, J(r’) x [V G(r,r’)] is perpendicular to the plane. Due to the “cross-product” with the 


normal vector, a zero value results for Lh[J(r)] = 2 n A x J/s’ J(r’) x [V G(r,r’)]dS’ . Hence, only the 
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2 n A x H lnc (r) term appears in the PO region. Our simulations for simple 2D geometries did yield 

this result, and hence, the hybrid scheme simplied, and did not require the use of an iterative 
technique for convergence. The detailed numerical values are, therefore, not included in this report. 

VIII. Results for the Two-Dimensional FDTD Scheme 

The two-dimensional FDTD code developed was applied for time-dependent simulations of 
scattering from biological cells. Since these cells function as dielectrics with finite conductivity, the 
application becomes more general than the “perfect-conductor” scenarios. There is also some interest 
in electromagnetic interactions with biological cells. Hence, some examples of light-scattering off 
biological cells using the 2D FDTD code were used & are presented here. 

As is well known, Mie theory has the two following fundamental limitations, (i) The incident 
energy has to be in the form of a plane wave. This has recently been circumvented by the 
development of a generalized Lorenz-Mie theory [44] that allows for the interaction between spheres 
and arbitrary incident beams, (ii) A spherical geometry and homogeneity of the scatter. While recent 
techniques have been developed to treat multi-layered spheres [45,46], the angular symmetry remains 
a constraint. The FDTD approach used here circumvents such restrictions, and allows for the 
inclusion of a heterogeneous structure, arbitrary shapes, and complex, non-concentric dielectric 
distributions. Having validated the FDTD method with Mie theory for the simple, homogeneous, 
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spherical case, results of the scattering patterns are now presented for the more realistic and complex 
cellular shapes. 

Different parts of a biological cell have different index of refraction as shown schematically in 
Fig. 6. Hence, cellular scattering patterns can potentially reveal various geometries of sub-cellular 
structures, their relative size and physical location relative to the center. For example, the angular 
scattering patterns can be obtained by changing the orientation of the incident optical light. If these 
were to be compared to the theoretical predictions for various shapes and sizes on the internal sub- 
structure, then the geometry, and size could be inferred. Similarly, changes in morphology (for 
example between healthy and malignant cells, or upon electromagnetic pulse exposure) can similarly 
be obtained. Some experimental work, in this regard, already exists. For example, Mourant et al. 
[47] showed that organelles smaller in size than the nucleus were responsible for significant 
scattering contributions. Measurements also indicate that the nucleus is responsible for low-angle 
scattering, while smaller organelles produced high-angle scattering. 



Fig.7. Cellular scattering patterns of cell with nucleus and without nucleus. 
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Fig. 7 shows FDTD results of the cellular scattering patterns for cell with nucleus and without 
a nucleus for an incident wavelength of 0.6328 pm, which correspond to a He-Ne laser source. In 
keeping with published reports [39], the refractive index of the nucleus and cellular cytoplasm were 
set to 1.39 and 1.37, respectively. The extra-cellular fluid was taken to be 1.35 [38]. The cell radius 
was 0.525 pm, while the nuclear radius was 0.263 |Jm. The patterns show that the nucleus brings 
about much more change at low angles than at higher angles. This is because the existence of a 
nucleus changes the forward scattering more significantly than backscattering. This prediction is in 
keeping with the experimental observations of Mourant et al. [47], Cells with mitochondria were 
also simulated. The scattering patterns obtained, for cells of radius 0.525 |jm are shown in Fig. 8. 
Both of the two cells had identical nucleii. As apparent from the curves of Fig. 8, the existence of 
mitochondria is predicted to lead to stronger backscattered intensity. Also, the scattering pattern is 
seen to change only at the higher angles. This trend also agrees well with measured data [47] 
suggesting the organelles can significantly change the higher-angle patterns. Thus, any 
morphological changes or electromagnetic field impact on the organelles (as in intra-cellular 
electroporation ), is likely to be diagnosed through such optical scattering. 



Fig. 8. Scattering patterns of cell with nucleus, with and without mitochondria. 
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All models used in literature implicitly assume that the nucleus is located at the center of a cell, and 
that its shape is spherical [48,49], However, this is a simplification, and in reality, the nucleus can 
have a number of different shapes and can be located anywhere inside a biological cell. The size, 
position and shape can all have an effect on the scattering patterns as demonstrated through the 
FDTD simulations. The results shown in Fig. 9 reveal the scattering patterns of cells with different 
nuclear sizes for a fixed outer cell size. The outer cell radius was again set at 0.525 pm, while the Ne- 
He laser characteristic wavelength of 0.6328 |Jm was used. With increasing size of the nucleus, the 
scattering intensity is predicted to increase as may be expected. Furthermore, the change in the 
scattering pattern is most pronounced at larger-angles. Hence, the backscattering configuration 
would be most sensitive for detecting such information. It is also apparent that as the size of the 
nucleus gets progressively smaller, there would not be much change in the scattering pattern and so 
the sensitivity would reduce. Hence, a size threshold is probably inherent from practical detector 
considerations. Furthermore, if sizes of the nucleus were the same, but a nucleus was located at 
different positions inside the cell, then again differences in the scattering pattern are predicted. This 



Fig. 9. Cellular scattering patterns of cells with different sizes of nucleus. 
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is shown through the results of Fig. 10. Two situations have been considered: one in which nuclei 
are positioned at the same radial distance from the center of the cell. However, their angular 
orientation relative to the incoming wave is different. A second situation is that of an identical 
nucleus positioned at the cell center. As seen from the results of Fig. 10, the scattering patterns for 
the former situation do not change significantly, and are all close to each other. However, they are all 
quite different from that for the center-located nuclear cell, and also have slightly higher intensities. 
The absence of sharp minima in the scattering pattern could, for example, be indicative of nuclei that 
were off-center in biological cells. Another case considered here is that of a nucleus located at 
different distances from the center of the cell. The results are shown in Fig. 1 1 with the relative 
orientation with respect to the incident wave being fixed. The scattering patterns once again reveals 
that changes are most likely to occur at higher angles (i.e., backscattering geometry). It is also 
apparent that the scattering pattern is affected by the relative distance of the nucleus from the center, 
but not, the absolute location. 



Fig. 10. Scattering patterns of cells with nucleus at the same distance from the centers of the cells. 
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IX. Results for the Three-Dimensional ADI-FDTD Scheme 

A three-dimensional (3D) ADI-FDTD simulation scheme was developed and implemented. 

The primary advantages of such an ADI-TDFD technique can be summarized as follows: 

(i) It allows for both a complete time-domain analysis and the steady-state result. The latter 
can result by running the simulations up to longer times. The transient behavior cannot 
be obtained from the Moment Method or PO-MoM schemes. 

(ii) The ADI-FDTD allows for any arbitrary geometry including sharp comers and complex 
features. 

(iii) It is significantly faster than the FDTD scheme since very large time steps can be used 
without the need to be restricted to the CFL condition. 

(iv) High frequency analysis for which the spatial discretization has to be small, can be carried 



Fig. 1 1 . Scattering patterns of cells with nucleus at different distances from the centers of cells. 
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out relatively easily as the time-step can remain large. Similarly, the technique is of 
advantage for large objects for which many grid points may be required. 

(v) Non-uniform, multi-grids can be employed for better resolution. Also, absorbing 
boundary conditions can be implemented for accurate analyses. 

(vi) Finally, dielectrics or perfect conductors or both can be studied with ADI-FDTD. 

The first set of results that are given below, emphasize the significant increase in time step that can be 
achieved through the ADI-FDTD scheme as compared to the traditional FDTD method. A simple 
electric field pulse was taken for simplicity. The traditional FDTD required a time step of about 
2.476x 10" ,7 s. Fig. 1 2a shows the electric field waveform ofan assumed source (single pulse), and 
also the time-dependent E-field at an observation point slightly “downstream”. The time step was 
chosen to equal the normal CFL (CFLN) value, and hence, was very small. Fig. 12(b) shows the 
numerical results of the E-field as a function of time, but with a time step chosen to be 1 000 times 
larger. The source E-field pulse was correspondingly chosen to be 1000 times longer. In Fig. 12(b), 
the E-field at the downstream point and the source waveform almost coincide since the time delay is 



Fig. 12a. Source E-field pulse and E-field at an observation point downstream obtained from an 
ADI-FDTD solution. The time step was chosen to equal the CFL condition. 
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Fig. 12b. Source E-field pulse and E-field at an observation point downstream obtained from an 
ADI-FDTD solution. The time step was chosen to equal 1000 times the CFL condition. 

very small compared on this new time-scale. The significant point, though, is that the field at the 
observation retains the same shape and does not cause any numerical instability despite the larger 
time step. Finally, a time step significantly larger ( 1 0 6 times) the CFL requirement was used to check 
for the numerical stability of the ADI-FDTD code. The time-dependent electric fields of the source 
pulse (correspondingly taken to be 10 6 times larger), and at the observation point are shown in Fig. 
12(c). Due to the large times, the delay between the E-fields of the source and at the observation 
point again cannot be discerned in the figure. However, both values were calculated and match the 
expected trend without any instability. This, conclusively demonstrates the successful 
implementation of an ADI-FDTD technique that can use time steps up to 1 0 6 times the CFL condition 
of traditional FDTD. 

A second example is now presented and discussed with the aim of demonstrating a successful 
comparison between ADI-FDTD and traditional FDTD. However, in this case absorbing boundary 
conditions are implemented through a perfectly matched layer proposed by Berenger [30]. A source 
from an incident plane is used and a thin dielectric wall placed between the source plane and the 
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Fig. 12c. Source E- field pulse and E-field at an observation point downstream obtained from an 
ADI-FDTD solution. The time step was chosen to equal 10 6 times the CFL condition. 


Observation pt. 



Fig. 1 3a. Geometry of a planar E-field source and an observation point separated by a thin dielectric. 
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observation point as shown in Fig. 13(a). For an incident sine wave, the electric field at the 
observation pointas obtained by the ADI-FDTD method in the absence of a dielectric wall is shown 
in Fig. 13(b). Due to the finite separation, there is a time lag between the two curves. Fig. 13(c) 
shows the E-field output at the observation point with ADI-FDTD and the traditional FDTD 
techniques. The time step for the ADI technique was 8 times larger. Though larger time steps could 
have been chosen from a stability standpoint, it would then not have been possible to plot the values 
over the 3.5 x 1 O' 10 s range of the graph. The close match between ADI-FDTD and the conventional 
FDTD is brought out clearly in Fig. 1 3(c). 

The use of multi-grid meshes for solving ADI-FDTD problems is brought out through a final 
example. Shown in Fig. 14 is the geometry under consideration. It consists of an incident source 
plane, a surrounding PML region for absorbing boundary conditions, four observation points, and an 
embedded cube of water with a finite conductivity (to facilitate absorption) and a dielectric constant. 
The water cube was assumed to have a dielectric constant of 1 .69, with a conductivity of 0.03 s/m. 
This cube was taken to have a material coating of conductivity 4.0 s/m and dielectric constant of 4.0. 



Fig. 13(b) Incident sine wave and the wave at the observation point without the dielectric wall. 
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Wavelength X = 0.6328jam. Cell size = 10.5*10.5*10.5 nm. 


conventional FDTD 

O AD I- FDTD 


lime (s) x1Q -io 

Fig. 13(c). Output wave at observation point calculated by conventional FDTD and ADI -FDTD. 
Time step for the conventional FDTD was 2.022 X 1 O' 12 s and that for ADI-FDTDwas 1.618X 10' u s. 
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Fig 14. (a) A water cube (o = 0.03s/m, er = 1 .69) coated with material (o = 4.0s/m, er = 4.0) 
(b) Multi-sized grid in FDTD calculation. 





The grid size was assigned as follows: Ax(i) = 1 .05x1 0' 7 m (in region I) ; Ax(i-l) x0.5 (In region 
II); lxl 0‘ 8 m (in region III) ; Ax(i-l) x2.0 (in region IV) ; and = 1.05xl0' 7 m (in region V). The 
spacings Ay(i) and Az(i) are defined in the same ways as Ax(i). Fig. 15(a) shows the output 
waveform at observation point 2, and compares the results of FDTD and ADI-FDTD. The time steps 
were 0.02022 fs and 2.022 fs, respectively. CPU times of 2985 s and 257 s were required for the 
FDTD and ADI-FDTD schemes, respectively. The savings with the ADI technique are clear. The 
output at observation points 1 and 2 are given in Fig. 1 5(b). 

X. Summary and Potential Future Work : 

A study into the problem of determining electromagnetic solutions at high frequencies for 
problems involving complex geometries, large sizes and multiple sources (e.g. antennas) has been 
initiated. Typical applications include the behavior of antennas (and radiators) installed on complex 
conducting structures (e.g. ships, aircrafts, etc..) with strong interactions between antennas, the 
radiation patterns, and electromagnetic signals is of great interest for electromagnetic compatibility 



Fig. 15(a). Output waveform at the close observation point. CPU times: 257 s for ADI-FDTD, and 
2985 s for conventional FDTD. 
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Fig. 15(b). Output at close observation point and far observation points 1 and 2. 


control. This includes the overall performance evaluation and control of all on-board radiating 
systems, electromagnetic interference, and personnel radiation hazards. 

Electromagnetic computational capability exists at NASA LaRC, and many of the codes 
developed are based on the Moment Method (MM). However, the MM is computationally intensive, 
and this places a limit on the size of objects and structures that can be modeled. Here, two approaches 
are proposed: (i) a current-based hybrid scheme that combines the MM with Physical optics, and (ii) 
an Alternating Direction Implicit-Finite Difference Time Domain (ADI-FDTD) method. The essence 
of a hybrid technique is to split the overall scattering surface(s) into two regions: (i) a MM zone 
(MMZ) which can be used over any part of the given geometry, but is most essential over irregular 
and “non-smooth” geometries, and (ii) a PO sub-region (POSR). Currents induced on the scattering 
and reflecting surfaces can then be computed in two ways depending on whether the region belonged 
to the MMZ or was part of the POSR. For the MMZ, the current calculations proceed in terms of 
basis functions with undetermined coefficients (as in the usual MM method), and the answer obtained 
by solving a system of linear equations. Over the POSR, conduction is obtained as a superposition of 
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two contributions: (i) currents due to the incident magnetic field, and (ii) currents produced by the 
mutual induction from conduction within the MMZ. This could effectively lead to a reduction in the 
size of linear equations from N to N - Npo with N being the total number of segments for the entire 
surface and Npo the number of segments over the POSR. The scheme would be appropriate for 
relatively large, flat surfaces, and at high frequencies. 

The ADI-FDTD scheme is relatively new in the literature and was first proposed capable in 
2000. It is capable of providing both transient and steady-state analyses. The primary advantage of 
this scheme is that large time-steps can be chosen without instability or inaccuracy, even though the 
spatial resolution is very small. Thus, computational speed enhancements of an order or more 
become possible. Hence, the ADI-FDTD method can be applied to problems involving large 
geometry, low wavelengths (i.e. high frequencies), and complex material characteristics such as 
composites and dielectric coatings. 

This report includes the problem definition, a detailed discussion of both the numerical 
techniques, and numerical implementations for simple surface geometries. Numerical solutions have 
been derived for a few simple situations. The advantage of the ADI-FDTD scheme has been 
demonstrated. Future work in this area might include: 

(i) Direct comparison of numerical results from the present MM-code at NASA LaRC with ADI- 
FDTD results for simple geometries. 

(ii) Calculations for complex geometries and/or composite material properties including dielectric 
coatings. 

(iii) Expansion of the ODU codes to include wires, wire-surface and wire-wire interfaces. 

(iv) Use of the codes in the choice of materials and design for electromagnetic applications. 
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